{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial - Sampling Problem\n",
    "---\n",
    "The final state of an annealing process is a Gibbs state, which is a resource we can use to produce some interesting probabilistic distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ising model\n",
    "\n",
    "A $3\\times 3$ spin lattice has a Hamiltonian like\n",
    "\n",
    "$$ H = \\sum_{\\mbox{i and j are neighbors}}Z_iZ_j .$$\n",
    "\n",
    "Under a thermal equilibrium, we define an order parameter \n",
    "\n",
    "$$ P = \\frac{1}{N}|\\sum_{i=1}^N(-1)^{x_i+y_i}S_i |.$$\n",
    "\n",
    "$S_i$ is the measured spin value on site $i$, $x_i, y_i$ are coords of this site, $N$ sites totally.  We will study how this parameter changes with inverse temperature."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set up the model\n",
    "\n",
    "We map this spin lattice to a 9-qubit annealer.  Each qubit plays as a spin on a site and is assigned a pair of coords. For example, the 3th qubit is located at $(1, 0)$, the 4th is at $(1, 1)$, and so on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from phalanx import Annealer\n",
    "\n",
    "beta = 1.0\n",
    "num_qubits = 9\n",
    "ism = Annealer(beta, num_qubits)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "connection = [\n",
    "    [0,1,0,1,0,0,0,0,0],\n",
    "    [0,0,1,0,1,0,0,0,0],\n",
    "    [0,0,0,0,0,1,0,0,0],\n",
    "    [0,0,0,0,1,0,1,0,0],\n",
    "    [0,0,0,0,0,1,0,1,0],\n",
    "    [0,0,0,0,0,0,0,0,1],\n",
    "    [0,0,0,0,0,0,0,1,0],\n",
    "    [0,0,0,0,0,0,0,0,1],\n",
    "    [0,0,0,0,0,0,0,0,0],\n",
    "]\n",
    "ism.set_connection(connection)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEuCAYAAADx63eqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABHlklEQVR4nO3de1yUZd4/8M89R2YAcTgICMMACqiAaaDmmnlCM7MUykNrKllP+2tPaes+e2zbfZ6nPbS7as/r9/u1W/srybZN2zxkmaVmlpYpKCmeRgWG4SgwAwgzMKfr9wcLD8jpnmGGe27m+369eInGzP3FcD5zXff3ui6OMcZACCGEBAiJ0AUQQgghI4mCjxBCSECh4COEEBJQKPgIIYQEFAo+QgghAYWCjxBCSECRCV0AIYSMRowx7Dtfhe1H9Whss8Fqd2KwxWMcB6jkUkQEK7AlJxW50+PAcdzIFRxAOFrHRwgh3lXTbMXm3cW4WNUMi83p9uNVcimmxodhx5ppiA1T+aDCwEbBRwghXnS23IT8nWfQYXfB4fL85VUm4aCUS1CQPxPZieFerJBQ8BFCiJecLTdhw+tnYLW7P8obiEouxa5NFH7eRM0thBDiBTXNVuTv9G7oAYDV7sTGnWdQ02z16vMGMgo+QggZJsYYnn2nGB12l0+ev8PuwubdxaAJOu+g4COEkGHad74KJdXNw7qnNxiHi+FCZTP2F1f55PkDDQUfIYQMA2MM24/qPeredIfV7sS2I3oa9XkBreMjhJBhKDSY0dhm4/31trpSmD97A7aaG2AOG2Rh4xCatRyhdz845GMb22woMpip0WWYaMRHCCHD8G6R0a2Gllvv/Rfay85DpomBOu1bsDdWwvTJK2g3XBjysVa7E3uKjMMpl4BGfIQQMiyF5eZBd2TpiTkdcN5uAABELHsWiqhE2BuNsNXegKO5bujHs87rkeGh4COEEA/ZnS4YTRbeX89JZQjNfgi3zx5A46GXIQ+Ph632JuTjkqBOnc3rOYxmC+xOF+RSmrDzFP3NEUKIhxpaOyBzM4DUKbMhDYuGreY62i4dByRSqFPuAafgtzWZTMKhobXDk3LJv1DwEUKIh+xOBnf2kXZaW3Dr3RfgbK5D9Lo/IH7zO1BEJ6H51D/QWnyY13NwHAe7kzo7h4OmOgkhZBAulws1NTUwGAwoLy+HwWDo/ry8rgnWRT8GJw/i9VyOpjowewcgkUEZmwpOJoc8QgtbzXXYG/g1rTDGIJfSqQ3DQcFHCAloDocDlZWV/QabwWCA0WiERqNBYmIidDodEhMTMXXqVDz00EOIT9DhkXcqeC9cl0doIQkKhav9Nure+QVkY2PQdvlzAIBSO4VfvS6GyBClx98voU2qCSGjXEdHByoqKgYMtpqaGkRHR/cKNp1O1/15QkICgoIGHtEt/PNnKG1o419P9TU0fb4Lttqb3ev4QqYtxZgZK3g9PjkyGJ/+aD7v65G+aMRHCBG1tra2PmHW8/PGxkbExcX1CraFCxd2fx4fHw+5XO7x9bMTNShrbOO9pEE5Pg3Ra//Lo2txXOf1yPBQ8JGAYHe60NDaAbuz8/5IZIiS2sFFoqmpadBga21tRUJCQq9gW758effnsbGxkEqlPqtvVZYWH1yo8fmWZUDnEUWrs7Q+v85oR8FHRiXGGAoNZrxbZERhuRlGkwUyqQQc17kI2OF0QRuuRnaiBquytMjWacC5055HvIIxhoaGhkGDzel09pqCTExMxKxZs7o/HzdunKD/77J1GkQEK2Cx+f7YoMgQJbJ0NOIbLrrHR0YVxhj2na/C9qN6NLbZYLU7B52C4rjOd9ERwQpsyUlF7vQ4CkAvcrlcqKur6w6x/oJNqVT2Cbae99g0Gv9/U7L3XCV+eaDEp6M+lVyK3+ZmIHd6vM+uESgo+MioUdNsxebdxbhY1ezRC5BKLsXU+DDsWDMNsWH8FhMHOofDgerq6gGDzWg0IiwsbMBg0+l0GDNmjNDfxrAxxrD2tdMoMph9cjSRTMIhS6fBO/92j9+/CRADCj4yKpwtNyF/5xl02F3DeuGRSTgo5RIU5M+kHfAB2Gw2GI3GAYOtpqYGUVFRAwZbQkIC1Gq10N/GiKhptiJn+wm0dXh/1BeslOLYlvmICeO3XpAMjoKPiN7ZchM2vH7GrR3yh6KSS7Fr0+gPP4vF0h1k/QVbQ0MDxo8fP2CwabVaKBQKob8Nv1FYbsJ6+ln0exR8RNR8/S776JZ5op72bGlp6fe+WtfnLS0tSEhIGDDYxo8fD5mMeuDcUVhuwkaaffBrFHxEtBhjWPPqaZyrCMz7KowxmEymXoF2Z7DZ7fZejSJ3Blt0dDQkElrW4W1d95svVDZ7NPpTK6TIjKP7zb5CwUdEa7R30jHGUFdXN2iwyeXyQYMtIiLCL0M7EDDGsL+4CtuO8OswBhikzInx4aF4bnEqVk6jDmNfoeAjosQYw31/PA6j2fdrp7QaFT7/8QKvvwg5nU5UV1cPGGwVFRUIDQ0dNNjCwsK8WhPxPsYYigxm7OlaU2q2QCbhwHEcGGNwuBi0GjUmRcjx3p+2wlB4fNAt0sjw0eQ9EaVCgxmNbbYhv67dcAF1//h5v/8tYtlmhEzNGfI5GttsKDKY3b7PYrPZUFlZOeA9tqqqKkRGRvYKtqysLDzyyCPdHZHBwcFuXZP4H47jkJ0Y3v3zM9guQpd3KnHo0CHk5eUJWfKoR8FHROndIiOveyfSMZEIzX64+/fM1o7WC58AAGSaWF7Xstqd2FNk7BN8VqsVFRUVAwZbXV1dd0dkV7DNnTsX69ev7+6IVCppl/1AI5dKBrxvt3HjRhQUFFDw+RhNdRJRcndH/C4thQdhPvpXKKInIPaJl3k/LlxmQ07HV72CrampCVqtdsCpyLi4OOqIJG65ffs2tFot9Ho9xo0bJ3Q5oxb9qySiY3e6YDRZ3H4cYwy3i94HAITyPAKmi9kuQ5gmHLnTp3cHW0xMDHVEEq8KDQ3FQw89hLfffhubN28WupxRi4KPiE5DawdkUgnsLve6Oa03zsBhroE0JBzBk+e69dgghQybvvsstZYTn8vPz8fWrVsp+HyI3q4S0bE7GTxpsLxdeAAAEDL9AXBS985f4zgOdifdFSC+t2DBAjQ2NuLChQtClzJqUfARUWGMoa6mGg6Hw63H2W6Vo91wAZxMgdDpyzy6rlxKa6qI70kkEqxfvx4FBQVClzJq0VQn8VsmkwkXL17s9VFSUoKQ0DAoHv8/gIT/4aIt/xrtBU+ZD6na/bVvDhdDZAh1YJKRsWHDBsybNw9/+MMfqEHKB+hvlAjOarXiypUrfUKutbUVGRkZyMjIQGZmJh577DFkZGQgMjLSra5Op6UZlssnAAChMx4e4qv7p9Wo6cR2MmLS0tKQlJSEjz/+GA8++KDQ5Yw6FHxkxDidTty8eRMlJSW9Aq6iogITJ05EZmYmMjIy8P3vfx+ZmZnQ6XQD7paSnahBWWPbEFtAdZKqw5Cwda/HdXNc5/UIGUlda/oo+LyP1vERr2OMoba2tntqsivgrly5gqioKGRmZnZ/ZGRkIC0tze2jbc6Wm7DxjTM+3aezi1ohxZtP0A75ZGSZzWYkJiairKwM4eH0s+dNNOIjw3L79u3ucOsZcoyx7nCbPXs2nn76aaSnp3vttO1snQYRwQpYbL7fqzMyRIksHY34yMjSaDS4//77sXv3bjzzzDNClzOq0IiP8GK323Ht2rU+jSZ1dXWYPHlyr1FcZmYmYmJifL6z/Gg/nYGQQ4cO4T/+4z9w+vRpoUsZVSj4SC+MMVRUVPRpNLlx4wYSEhK6pye7Am7ChAmQSvl3V3q71rWvnUaRITDP4yOjn8PhgFarxWeffYa0tDShyxk1KPgCWGNjY58pypKSEoSEhPS6B5eZmYkpU6ZApfK/XUt8fQL7sS3zERNGR8QQ4WzduhUKhQK//e1vhS5l1KDgCwBWqxWXL1/uE3JdywXuDLmIiAihS3ZLYbkJ618/49FJ1wNRyaXYtYkaWojwLl68iGXLlqG8vFyw2ZXRhoJvFOlaLtBz9HbncoGeITfYcgGxKSw3YePOM+iwu4Y17SmTcFDKJSjIp9Aj/uPuu+/GSy+9hJycoc+PJEOj4BOhnssFeoZcf8sFMjMzkZqa6vZyATGqabZi8+5iXKhs9mj0p1ZIkRkXhh1rptFm1MSvvPzyyygsLMSuXbuELmVUGJHgG+zEYTK4nssFeoZcz+UCXR/eXC4gVowx7C+uwrYjejS22WC1Owdf5M5cUCvliAhW4LnFqVg5LW7UjILJ6FFfX4+UlBRUVFQE/L9xb/BJ8DHGUGgw490iIwrLzTCaLJBJJeA4gDHA4XRBG65GdqIGq7K0yNZpAv7FxmazQa/X9+mmvHXrFqZMmdKnm3IklguIGWMMRQYz9nT9DJotkEk4cBwHxhgcLoboYBkqio7jn3/cihmJ4fT3SfzaihUrsGLFCmzatEnoUkTPq8HHGMO+81XYfpTfu22O62wiiAhWYEtOKnKnj/5324wxGAyGPo0mPZcL9Aw5IZcLjCYDzTpMnjwZBQUFmDlzptAlEjKoffv2YceOHThx4oTQpYie14Kv6/7KxapmjxYUq+RSTI0fXfdXupYL9Ay5O5cLdH1MnjzZL5cLjHY///nP4XK58Pvf/17oUggZlM1mQ1xcHM6cOYOkpCShyxE1rwTf2XIT8gO4o67ncoGeQXfncoGukZzYlguMZoWFhVi3bh2uXr066mcbiPj94Ac/QGRkJF544QWhSxG1YQff2XITNgTIGqo7lwt0fRiNRqSkpPQJuYSEBHox9XOMMeh0Ohw+fBhTpkwRuhxCBlVYWIg1a9bgxo0b9NoyDMMKPl/vmnF0yzxBpj0ZY6ipqenTTXnlyhVER0f3aTQJlOUCo9UPf/hDREdH4xe/+IXQpRAyKMYYMjIy8Je//AVz584VuhzR8jj4GGNY8+ppnKsQ9z6JLS0tKCkp6RNyHMf1aTSh5QKj0/Hjx7F161YUFRUJXQohQ3rppZeg1+vxt7/9TehSRMvj4BPbzvg2m637dIGeIVdfX9/rdIGukKPlAoHD4XAgNjYWhYWF0Ol0QpdDyKCqq6uRnp6OqqoqqNVqocsRJY+CjzGG+/54HEaz789C02pU+PzHC3iHUM/lAj1D7s7lAl0fycnJtFyAYNOmTbjrrrvw7LPPCl0KIUNaunQp1q9fj3Xr1gldiih5FHzunn5tufYlmr96F/YGAyCVQRGViKhHfwVpUMiQjx3s9OueywW6Pi5duoTQ0NA+jSa0XIAM5uDBg/jzn/+Mzz77TOhSCBnSP/7xD7zxxhv45JNPhC5FlDwKvn9/7xu8W1Q5+FZQ/9J2+QQa3v8jIJVDnXoPJHIVOmr0GLf615CFRg5dIAfk3hWLdRNZn5CzWCzdAdfzV1ouQNzV3t6OmJgYXL9+HVFRUUKXQ8igrFYr4uLicOHCBcTH0yHJ7vIo+Bb++TOUNrQN+XWMMVS9sgnOlnpEP/ZbBOmmelSkw1SF8K/+d5/7cLRcgHjTqlWrsHTpUjz55JNCl0LIkJ5++mkkJyfjpz/9qdCliI7bwWd3ujDlV4dh59HJaTdVofrV74CTKaFMyECH8RKkwRqMmbECoVnLeV9TJuFw5T+W0sbWxKf+8Y9/4O9//zs++OADoUshZEinTp3CU089hcuXL9MAwE1uJ0lDawdkPAPIaWkBADBHBxxNdVBPuhfO1kaYjvwFFv1XvK8pl3JoaO1wt1RC3LJs2TJ8/vnnuH37ttClEDKkb33rW3A4HDh79qzQpYiO28FndzLwfXMhVf/PmrfIh55D5IObETx1MQDAcv1r3tfkOA52Jx0bSHwrLCwMc+bMwUcffSR0KYQMieM4bNiwAQUFBUKXIjpuB59cyvFqagEAWdg4cMr+15lIFPw7LBnr3FGfEF/Lzc3F3r17hS6DEF7Wr1+P3bt3o6ODZsTc4XbwRYYo4XC6eH0tJ5VjTPYKAEDDB9vQ8OEOtF04AnASBKfP531Nh4shMkTpbqmEuG3FihU4fPgwvZAQUUhMTERmZibdl3aTByM+CbTh/HcLCJuzFmPueRSsvQ2Wq19AHqnDuEefh3J8Gu/n0GrU1NhCRkTXXqzHjh0TuhRCeNm4cSN27twpdBmi4vN1fMPFccCqrHi89Mhdvr8YIQC2bduGK1eu4LXXXhO6FEKG1Nraivj4eFy7dg3R0dFClyMKHg2jVmVpoZKPzDZfKrkUq7O0I3ItQoDO+3wHDhyA0+m7fWgJ8ZaQkBCsWLECb7/9ttCliIZHwZet0yAieGSO4YkMUSJLpxmRaxECAElJSYiLi8OpU6eELoUQXjZu3EjdnW7wKPg4jsOWnFSoFb4d9ankUmzJSaHFmWTE5ebmYt++fUKXQQgv8+fPh9lsxjfffCN0KaLgccdI7vQ4ZMaFQSbxTSjJJBymxodh5bQ4nzw/IYPpCr5hnNNMyIiRSCRYv349jfp48jj4OI7DjjXToJT7pttSKZfg5TXTabRHBJGRkQGZTIbz588LXQohvGzYsAFvv/027Ha70KX4vWGlVmyYCgX5M73e6KKSS1GQPxMxYUFefV5C+OI4jqY7iaikpqYiOTkZH3/8sdCl+L1hD9eyE8Oxa9NMBCulw572lEk4BCul2LWp//P3CBlJeXl5FHxEVGhNHz8erePrT02zFZt3F+NCZTOsdvfbwNUKKTLjwrBjzTTEhtGBsUR4LpcL8fHx+Oyzz5Camip0OYQMqampCTqdDmVlZQgPp8HDQLx2gy42TIV3/u0e/DY3A1qNCmqFdMjNrDmuM/C0GhVeXJmBd/7tHgo94jckEglWrFhBoz4iGmPHjsUDDzyAd955R+hS/JrXRnw9McZQZDBjT5ERheVmGM0WyCQcOI4DYwwOF4PdXIvF05PxnZypyNJpqImF+KVPPvkEv/rVr3D69GmhSyGEl48++gi//vWv8fXX/E/ACTQ+Cb472Z0uNLR2wO7sPGUhMkSJ7z3zvzBlyhRs3rzZ15cnxGM2mw0xMTG4ePEi4uJoaQ3xfw6HAwkJCfj0008xadIkocvxSyOy87NcKkFsmAoJ4WrEhqkgl0qQk5ODo0ePjsTlCfGYQqHAgw8+iAMHDghdCiG8yGQyrFu3jtb0DUKwIw8WLlyIL774gtacEL9HyxqI2GzcuBG7du2i/WYHIFjwRUZGYsKECThz5oxQJRDCy/3334+vv/4aJpNJ6FII4SUjIwPR0dH49NNPhS7FLwl6yN2iRYtoupP4veDgYCxatIgO+ySiQhtXD0zQ4MvJyaEDP4ko0HQnEZvHHnsMBw8eREtLi9Cl+J0R6eocSFtbG6Kjo1FbW4uQkBChyiBkSCaTCUlJSaiurkZwcLDQ5RDCS25uLpYvX44nn3xS6FL8iqAjvuDgYGRnZ+Pzzz8XsgxChhQeHo4ZM2bQPohEVGi6s3+CBh8AWtZARIP27iRis2zZMly9ehWlpaVCl+JX/CL46D4fEYMVK1bgww8/pCU4RDQUCgXWrl2LN998U+hS/IrgwZednQ2DwYC6ujqhSyFkUHFxcUhNTcVnn30mdCmE8LZx40a8+eabcLlcQpfiNwQPPplMhnnz5tF6EyIKubm52Lt3r9BlEMLb3XffDbVajZMnTwpdit8QPPgAmu4k4pGbm4sDBw7Qu2ciGhzHUZPLHfwm+I4cOQIBV1YQwktqairCw8Np53siKo8//jj27t2LtrY2oUvxC34RfJMmTYLdbsfNmzeFLoWQIdFidiI2sbGxmD17Nv3c/otfBB/HcTTdSUSj6z4fzVAQMaHpzv/hF8EH0Ho+Ih7Tp0+Hw+FASUmJ0KUQwtuKFStw7tw5GI1GoUsRnN8E36JFi/Dpp59S0wDxexzH0XQnEZ2goCA8+uijeOutt4QuRXB+E3xxcXEYN24ciouLhS6FkCFR8BEx6pruDPRper8JPoCmO4l4zJkzB1VVVSgrKxO6FEJ4mz17NlwuV8Cfg+pXwUfn8xGxkEqlePjhh2nUR0SF4zhs2LAh4JtcBD2W6E5NTU3QarWor69HUFCQ0OUQMqhDhw7hd7/7Hb744guhSyGEN4PBgKysLFRWVgbs66xfjfjGjh2L9PR0fPXVV0KXQsiQFi1ahIsXL9I+s0RUdDod7rrrLhw8eFDoUgTjV8EH0HQnEQ+lUomlS5fi/fffF7oUQtwS6Gv6/C74qMGFiAltWk3EKC8vD6dOnQrY2Qq/uscHAO3t7YiKioLRaMTYsWOFLoeQQd2+fRtxcXEwGo0ICwsTuhxCeMvPz8fUqVPx3HPPCV3KiPO7EV9QUBC+9a1v0ZlnRBRCQ0Nx33334dChQ0KXQohbAnm60++CD6D7fERcaDE7EaN58+ahubm516YhdqcLNc1WVJgsqGm2wu4cnTtp+d1UJwCcO3cO69atw5UrV4QuhZAh1dfXY+LEiaitrYVKpRK6HEJ4++Uvn4fBIkXsnFwUlpthNFkgk0rAcQBjgMPpgjZcjexEDVZlaZGt04DjOKHLHja/DD6Xy9W9fVl8fLzQ5RAypHnz5mHr1q146KGHhC6FkCExxrDvfBVe+ugSasytkCiCMFgScBygkksREazAlpxU5E6PE3UA+uVUp0QiwYIFC+iYIiIaeXl5NN1JRKGm2Yq1r53GLw+UoLbVAU4+eOgBnaM/i80Jo9mKX+wvwdrXTqOm2ToyBfuAXwYfQMsaiLisXLkSBw8ehMPhELoUQgZ0ttyEnO0nUGQww2JzevQcVrsTRQYzcrafQGG5ycsVjgy/Dr5jx44F/C7iRBx0Oh0SEhJo+zLit86Wm7Dh9TNo63DC4Rre66rDxdDW4cT618+IMvz8NviSk5OhUCiowYWIBnV3En9V02xF/s4zsNo9G+UNxGp3YuPOM6Kb9vTb4OM4jqY7iajk5eVh//79NEtB/ApjDM++U4wOu2+WJnTYXdi8u1hUP/d+G3zA/0x3EiIGkydPhkqlQlFRkdClENJt3/kqlFQ3D3t6cyAOF8OFymbsL67yyfP7gl8H38KFC3HixAlqGCCiwHEcTXcSv8IYw/ajeo8bWfiy2p3YdkQvmlGfTOgCBjNu3DjodDqcPXsWs2fPFrocQoaUm5uL/Px8vPjii0KXQggKDWY0ttl4f33t33+KDmNJrz+TRyZg/FP/d8jHNrbZUGQwIzsx3O06R5pfj/gAmu4k4jJjxgzcvn0bV69eFboUQvBukdGjhpbQ7Ie7P4IzFvJ6jNXuxJ4io9vXEoIogo8aXIhYSCQSrFy5kqY7iV8oLDcPuTi9P+E5T3d/hN3zKK/HMNZ5PTHw++CbO3cuCgsL0dbWJnQphPBC9/mIP7A7XTCaLB491rh9DSq2r0HdP36Ojho9/8eZLaLY2Nrvgy8kJAR33303Tp48KXQphPBy33334ebNmzAaxTHtQ0anhtYOyKTuvcRLFCqoJsyAevJcyMZEod1wAbd2/wrOVn4jOZmEQ0Nrhyfljii/Dz6ApjuJuMjlcixfvhz79+8XuhQSwOxOBnf3kY569FcYt+oFRCz9PmLzt0M6Zhxc7a1or7jA6/Ecx8Hu9P/OTlEEH53PR8SGNq0mQpNJOk+64ctlb4ezdYDtx3gmKGMMcqn/n9rgl8cS3clutyMyMhI3btxAVFSU0OUQMiSr1YqYmBjcvHkTkZGRQpdDRrHW1lbo9Xpcu3YN165d6/5cf+Mmwr9TAE7Kb9Wao6kOVa99B0G6uyAbE4WOqquw15dDEjwW45/8P5Cqw4Z8DrmUw+XfLIXczSnWkebX6/i6yOVy3HfffTh+/DhWr14tdDmEDEmlUiEnJwcHDx7EE088IXQ5ROScTifKy8v7hNu1a9dgNpuRkpKC1NRUpKWl4YEHHsDmzZuRmpqKvP9XjNIGfo2BElUoQjIWot1wAR0VF8Ep1VCl3IOx963nFXoAoNWo/T70AJEEH/A/050UfEQscnNzsWfPHgo+wltDQ0O/4VZWVobo6OjucEtPT0deXh5SU1Oh1WohkfQfNtmJGpQ1tvFa0iBRqhHxwA89rp3jOq8nBqKY6gSAkpISPPzwwygtLRW6FEJ4aWpqQkJCAqqrqxESEiJ0OcRPtLe348aNG72CrSvonE4n0tLSuj+6gi4lJQUqlcrta50tN2HjG2d8vmUZAKgVUrz5xExR7NwimhFfeno6LBYLSktLkZycLHQ5hAxp7NixmD17Ng4fPoxHH+W3CJiMDowxVFZW9htu1dXVSExM7A63e++9F08++STS0tIQFRUFzt1WzEFk6zSICFbAYvP9sUGRIUpk6cQx4hNN8HUdU3Ts2DEKPiIaubm52Lt3LwXfKNXS0tJvuF2/fh2hoaG9Rm6LFy9GWloakpKSIJONzEsvx3HYkpOKXx4o8emoTyWXYktOildD25dEM9UJAG+88QYOHz6M3bt3C10KIbzU1NRgypQpqKurg0KhELoc4gGHw4GysrJ+7721tLR0T0d2/dr1+ZgxY4QuHUDn6HPta6dRZDD75GgimYRDlk6Dd/7tHgo+X6ioqEBWVhbq6uoGvJlLiL+ZM2cOnn/+eSxdulToUsgAGGOor6/vN9zKy8sxfvz4XsHWFW5xcXGieC2qabYiZ/sJtHV4f9QXrJTi2Jb5iAkL8vpz+4qogg8A0tLSsHv3bkybNk3oUgjh5U9/+hOuX7+Ov/71r0KXEvCsViuuX7/eJ9z0ej04jusVbF3hNnHiRAQFiedFfSCF5Sasf/2MR6c1DEQll2LXJnE0tPQkuuD77ne/i+TkZGzdulXoUgjh5ebNm/jWt76F6upqSKVS2J0uNLR2wO7s3OUiMkQpirVPYuFyuWA0GvsNt9raWiQnJ/fbORkIGw0UlpuwcecZdNhdw5r2lEk4KOUSFOSLL/QAEQbf3r178eqrr+Lw4cNCl0IIL4wxpM9/GDNW/wDGdgWMJgtkUgk4rvMoF4fTBW24GtmJGqzK0iJbpxHNvRIhNTU1DdhYotFo+gRbWloadDrdiDWW+KuaZis27y7Ghcpmj0Z/aoUUmXFh2LFmGmLD3F9i4Q9EF3xmsxk6nQ719fVQKpVCl0PIgBhj2He+CtuP6lFjboWDcQA38MiO4zqnjiKCFdiSk4rc6XEBH4B2ux2lpaV9wu3atWuwWCwDNpbQusnBMcawv7gK247o0dhmg9XuHHSROwdApej82XxucSpWThP3z6bogg/oPOX6T3/6E+bNmyd0KYT0q+td9cWqZo/ayFVyKabGi/tdNV+MMdTV1fUbbhUVFYiPj+8TbmlpaYiNjRX1i68/YIyhyGDGniIjCsvNMJotkEk4cBwHxhgcLga1ywKN04w/PfMIskbJbIQog+9nP/sZZDIZ/vM//1PoUgjp42y5Cfl0H6WPtra2ARtLFApFv+GWnJxMMzsjqL/7z0Vnz+A73/kOvvnmG6HL8xpRBt+xY8fw/PPP48svvxS6FEJ6OVtuwoYA7pxzOp2oqKjoE27Xrl1DQ0MDJkyY0G9jSXi4/39vgcrhcCAqKgqXL19GbGys0OV4hSiDr729HVFRUaisrERYGL9dwwnxNV+vlTq6ZZ7fTHuaTKZ+R243b95EREREv+GWkJAAqVQqdOnEA4888ghWrlyJ9evXC12KV4iyvSkoKAizZs3CiRMn8PDDDwtdDiFgjOHZd4rRYed/8Kc7OuwubN5dPKK7Y3R0dODmzZv9nvXW0dHRK9xWr17dvZlycHDwiNRHRs6SJUtw5MiRURN8ohzxAcDvf/971NTU4OWXXxa6FEKw91zliOyH+NvcDOROj/faczLGUFNT029jSWVlJRISEvpdFhAdHT0qmhwIP2VlZZg9ezZqampGxf93UY74ACAnJwcbN24UugxCwBjD9qN6nx/9YrU7se2I3qNW8gFP6dbroVarewXbggULujdTpv1FCQAkJSUhJCQEFy9exNSpU4UuZ9hEG3zTp09HTU0NqqurMX78eKHLIQGs0GBGY5vNrce0XT6Bhvf/CAAIzX4Y4TlP83pcY5sNRQZzv40unp7SPXbsWLdqJ4FpyZIl+OSTTyj4hCSVSrFgwQIcO3Zs1Mw7E3F6t8joVheno6UBpo//LyCRAi73RolWuxNvnrqBjip7n67J/k7pzs3NRVpa2qCndBPCx5IlS/DKK6+Miu0iRXuPDwBeeeUVfP3119i5c6fQpZAAtvDPn6G0oY3X1zLGcOudX8DZ1gR5lA6WK1+4NeIDAKe5GuPP/81rp3QTwkdzczPi4+Nx69Yt0f+ciXbEBwCLFi3Ciy++CMbYqLjhSsTH7nTBaLLw/vrbZw+gvfIyYjdsQ8vZAx5dMygiDie//Io2tiYjKiwsDFOnTsXJkyexePFiocsZFlH/y0lJSYFEIsG1a9eELoUEqIbWDsh4BpCtvhzmEwUYO/dxKKKTPb6mTMqhobXD48cT4qmuZQ1iJ+rg4zgOOTk5OHbsmNClkABldzLwnWywXPsScDrQXnERt979DdoNnVtAWa9/DfNnO3lfk+M42J2ivUNBRKyrwUXsRD3VCXROd/7zn//E9773PaFLIaOczWZDWVkZrl+/Dr1eD71ej8vl1WjLzAcn47GfJGMAGNpLi3r9saO5Dh1VV3nXwVjnPoqEjLQZM2bAYDCgtrYWMTExQpfjMVE3twBAbW0tJk+ejPr6+oA/Z4sMn8vlQmVlZXew9Qw5o9GI+Ph4pKamdi8NSJ6Ygh98bofDgw1bGj7YjraSY243t8ilHC7/Zind4yOCyMvLQ15eHh5//HGhS/GY6JMiJiYG8fHxOHfuHGbOnCl0OUQEGGNoaGjoDrSeIXfjxg1oNBqkpqZ2fyxYsKAz5JKT+13QnVDCv6vTG7QaNYUeEUzXfT4xB5/oR3wAsGXLFkRFReHnP/+50KUQP9LS0oLr16/3GrV1BRzHcb3Cretj4sSJbh9i+u/vfYN3iyoHPcjTWzgOWJUVj5ceucv3FyOkH6WlpZgzZw6qq6tF200/KoLvgw8+wLZt2/Dpp58KXQoZYT03Ur4z4FpaWrqnJLt+7fqIiIjwWg1ny03Y+MYZn29ZBgBqhRRvPiGOI4rI6DVhwgTs378fmZmZQpfiEdFPdQLAvHnzsHbtWlgsFqjVaqHLIV7WdcbbnaM2vV6P6upq6HS67kDLzs7GY489htTUVIwfP35EdivJ1mkQEayAxWb1+bUiQ5TI0ml8fh1CBtPV3UnBJ6DQ0FBMmzYNp06dEv3CykDFGENtbW2fUZter0dZWRnGjRvXa9T2wAMPIDU1FTqdDnK5XNDaOY7DlpzUETmdYUtOiminl8josWTJEvz1r3/Fj370I6FL8ciomOoEgBdeeAHt7e34wx/+IHQpZBBNTU39dkzq9XqoVKruYOsZchMmTPD7kTxjDGtfO40igxkOl/f/SckkHLJ0mhE9j4+QgTQ1NUGr1aK+vh5BQUFCl+O2UTHiAzqPKdq8eTOAzm2kGlo7YHd2rneKDFFSF9wIslqtuHHjRr9Tk1artde9tuXLlyMlJQUpKSnQaMQ7hcdxHHasmYac7Sfg8MEJ7Eq5BC+vmU6hR/zC2LFjkZmZiZMnTyInJ0foctw2KkZ8jDF8deMW8ra+hLQ5y1Dd3LmNFMd1rhl2OF3QhquRnajBqiwtsnUaegEZJofDgfLy8j7Tknq9HvX19UhKSurTMZmSkoKYmJhR/XdfWG7C+tfPuHVaw1BUcil2baKGFuJffv3rX8NiseCll14SuhS3iTr4GGPYd74K24/q0dhmg6XDDnADj+w4rvNFJCJYgS05qcid7v6BnoHE5XKhurq636lJg8GA8ePH99sxmZCQAKlUKnT5giksN2HjzjPosLuGNe0pk3BQyiUoyKfQI/7nyy+/xPe+9z2cP39e6FLcJtrgq2m2YvPuYlysavaooUAll2JqfBh2rJmG2DBxH7ExXI2Njf1OS16/fh1jxozpM2rrWswtxrn9kdL183mhstmj0Z9aIUVmHP18Ev/lcDgQFRWFq1evIjo6Wuhy3CLK4DtbbkI+vaN2S2tra5/7bl3h5nQ6B1zMPWbMGKFLFy3GGPYXV2Hbkc4ZCavdOfgid5cL6iA5IoIVeG5xKlZOoxkJ4t9yc3Px6KOPYt26dUKX4hbRBd/ZchM20D2UftlsNpSWlvbbMWk2mzFx4sR+pyYjIyPpBdaHGGMoMpixp8iIwnIzjGYLZBIOHMeBMQaHi0GrUaPs7DH816YH8e3Fs+j/BxGFV155BadPn0ZBQYHQpbhFVMFX02xFzvYTaPNB11ywUoqjW+b5/bSSy+WC0Wjsd2qysrISWq2236nJ+Pj4EVnMTYY2UNfxT37yE8hkMrz44otCl0gILzdv3sTcuXNRVVUlqjdrogk+xhjWvHoa5ypG/zopxhhu3brV72Lu0tJSRERE9Bm1paamIjExsd9NlIk4nDlzBhs2bMCVK1cE/xkkhK/k5GS8//77yMjIELoU3kSzjm/f+SqUVDf7JPQAwOFiuFDZjP3FVcidHu+Ta9ypubm5VyNJz4CTy+W9Rm1d23BNnDgRwcHBI1IfGVkzZsxAW1sbrly5gilTpghdDiG8dG1fJqbgE8WIjzGG+/54HEaz7/dC1GpU+PzHC7z2jru9vb17E+U7pyZbW1v7jNy6FnN7cxNlIh7PPvssIiMj8fzzzwtdCiG87N27F6+99ho++ugjoUvhTRTB5+7u9w0fbEN7eTGc1hZIFGooYiZCM28jFDEThnysJ7vfOxwOGAyGfqcma2trkZiY2G/XZGxsLE1pkV4+//xzPPvss6JcG0UCU1NTExISEnDr1i3RLHESxVTnu0VGt7o4Hc23oEzIhESpRrvhAtrLzuFWoxHx331jyMda7U7sKTL2CT7GGGpqavqdliwrK0NMTEyvZpIHH3ywexNlOhme8NV1zllpaSmSk5OFLoeQIY0dOxbp6ek4deoUFi1aJHQ5vIjiFbmw3OzWIZ8x637f/XlH7Q3U7twM5+1GMKcDnHTwb5kx4OTVauzadaHPYm61Wt1rxDZ79uzuTZRVKv/uBiXiIJVKsXLlSuzduxdbt24VuhxCeOm6zyeW4PP7qU6704UpvzoMu5tNLS1FB2FvMKLd8A0cpiqMmZUHzYJN/B7scuIe425MSk3pde8tLCzMg++AEPd8/PHH+PWvf42vvvpK6FII4eXUqVP4wQ9+gHPnzgldCi9+P+JraO3ccNrucm/tnuXqKXQYSwAA0tBIKOP4d8mplHJsf+Vvfr+mj4xOCxYswLe//W1UVlYiPn5kOowJGY6ZM2eitLQUt27dwrhx44QuZ0h+v6LZ7mTwpP8jZt3vkbB1L6Lyfglnqwn1+38HR/MtXo/lOA52p18PhMkoplAosHz5cuzfv1/oUgjhRS6XY/78+Th69KjQpfDi98Enl3Ju3d9z2TvA/jU65GQKqJKzwCmCAJcTjqZaXs/BWOeOGoQIJS8vD++9957QZRDCW9d9PjHw+6nOyBAlHE4X76+3VV9Dw8E/QalNhyQoBB3GS2AdFkjUYVBED72cAehczB4ZovS0ZEKGbcmSJdiwYQPq6+sRFRUldDmEDGnJkiV48cUXwRjz+2VaIhjxSaANV/P+emloBGSa8WgvK0brN0fgam+FetK9iH7sRUiC+O14otWo6cR2IiiVSoX7778fBw4cELoUQniZMGEClEolLl++LHQpQ/L7ER8AZCdqUNbYxmvKUx4e12s5g7s4rvN6hAjtkUceQUFBAZ566imhSyFkSBzHYfHixfjkk0+Qnp4udDmDEsWwZlWWFir5yJzorZJLsTpLOyLXImQwy5Ytw8mTJ9HU1CR0KYTwIpb7fKIIvmydBhHBI3PqQGSIElk6GvER4YWGhmL+/Pn44IMPhC6FEF4WLlyIU6dOoaOjQ+hSBiWK4OM4DltyUqFW+HbUp5JLsSUnxe9vzJLA8cgjj2Dv3r1Cl0EILxqNpnv7Mn8miuADgNzpcciMC4NM4ptQkkk4TI0Pw8ppcT55fkI88dBDD+HYsWNoa2sTuhRCeOm6z+fPRBN8HMdhx5ppUMp9U7JSLsHLa6bTaI/4lfDwcMyaNQuHDx8WuhRCeBHDfT7RBB8AxIapUJA/0+uNLiq5FAX5MxETJo4jNUhgeeSRR2gxOxGNWbNmobS0FPX19UKXMiBRBR8AZCeGY9emmQhWSoc97SmTcAhWSrFrk3vn7xEyklasWIGPPvrI7xsGCAE6ty+bN2+eX29fJrrgAzrD7+iWecjSaTwe/akVUmTpNDi6ZR6FHvFrMTExyMzM9OsXEkJ68vfpTr8/lmgwjDHsL67CtiN6NLbZYLU7B13kzoEBTjviI8Pw3OJUrJwWR/f0iCjs2LEDFy5cwOuvvy50KYQM6fr161iwYAGMRqNfvsaKOvi6MMZQZDBjT5ERheVmGM0WyCQcOI4DYwwOF4NWo8Z0bRje+o/v4qv3/06nWxNRqaiowN13343a2lrIZKLYcIkEMMYYkpOT8eGHH2LKFP5Hwo2UUfEviOM4ZCeGd09Z2p0uNLR2wO7sPGUhMkTZvfem49RMFBQU4De/+Y2QJRPiloSEBCQlJeHEiROiOeWaBK6e25f5Y/CJ8h7fUORSCWLDVEgIVyM2TNVrw+knnngCO3fuhMvF/8QHQvwBLWYnYuLP9/lGZfANZtq0aQgPD8enn34qdCmEuCUvLw/79u2jN21EFBYuXIiTJ0/6ZTdywAUfAGzatAlvvPGG0GUQ4pbU1FRERETgq6++EroUQoYUHh6OKVOm4MsvvxS6lD4CMvi+/e1v48MPP6Rd74no0HQnERN/3b4sIIMvIiICS5YswTvvvCN0KYS4JS8vD++99x5GQTM2CQD+ep8vIIMP6GxyoelOIjaZmZmQy+U4f/680KUQMqR77rkHN27c8LvtywI2+JYsWYLKykpcunRJ6FII4Y3jONq7k4iGXC7H/PnzcezYMaFL6SVgg08qlWLjxo006iOik5eXR/f5iGj4432+UbFzi6f0ej3mzp2LyspKyOVyocshhBfGGBISEvDxxx/75eJgQnrS6/VYuHChX21fFrAjPqCzPTwlJQWHDh0SuhRCeOM4rrvJhRB/l5KSAplMhqtXrwpdSreADj6A1vQRcaJlDUQsOI7zu+7OgA++VatW4bPPPkNdXZ3QpRDC25w5c1BdXY3S0lKhSyFkSP52ny/ggy80NBS5ubl46623hC6FEN6kUilWrlxJoz4iCosWLcIXX3zhN9uXBXzwAZ1r+l5//XVaFExEhZY1ELEIDw/H5MmT/Wa7PQo+AHPnzkVHRwfOnj0rdCmE8DZ//nzo9XpUVVUJXQohQ/Kn6U4KPnTefKWdXIjYKBQKLF++HPv27RO6FEKG5E8NLgG9jq8no9GIu+66C1VVVVCpVEKXQwgvBw4cwI4dO3D8+HGhSyFkUDabDVFRUbh58yYiIyMFrYVGfP+i1Woxc+ZMevdMRGXJkiU4f/683+2FSMidFAoF5s2b5xfbl1Hw9dDV5EKIWKhUKtx///04cOCA0KUQMiR/uc9HwdfDihUrUFxcDIPBIHQphPBGe3cSsei6zyf0HTYKvh6CgoKwdu1aFBQUCF0KIbwtW7YMJ0+epIOVid9LTU2FRCLBtWvXBK2Dgu8OXd2dLpdL6FII4SU0NBTz58/HBx98IHQphAzKX7Yvo+C7w913340xY8bgxIkTQpdCCG+0dycRC3+4z0fLGfqxY8cOFBUVYdeuXUKXQggvJpMJSUlJqK6uRnBwsNDlEDKgxsZGJCUloaGhAQqFQpAaaMTXj3Xr1uHgwYNobm4WuhRCeAkPD8esWbNw+PBhoUshZFARERGYNGmSoNuXUfD1IyoqCosWLcKePXuELoUQ3mjvTiIWQt/no+AbAK3pI2KzYsUKfPTRR36zAz4hAxH6Ph8F3wCWLl0Kg8GAK1euCF0KIbzExMQgMzMTR48eFboUQgY1e/ZsXLt2DY2NjQAAu9OFmmYrKkwW1DRbYXf6tquemlsG8ZOf/ASMMbz00ktCl0IILzt27MCFCxdotoL4NcYYFq19GhHZD6JRMhZGkwUyqQQcBzAGOJwuaMPVyE7UYFWWFtk6DTiO89r1KfgGcfXqVSxYsABGoxEymUzocggZUkVFBe6++27U1tbSzyzxO4wx7Dtfhe1H9ahtaoPdBYAbeOKR4wCVXIqIYAW25KQid3qcVwKQpjoHMWnSJCQlJVGnHBGNhIQEJCcn0zpU4ndqmq1Y+9pp/PJACYxmK+xMMmjoAZ2jP4vNCaPZil/sL8Ha106jptk67Foo+IZATS5EbGjvTuJvzpabkLP9BIoMZlhsTo+ew2p3oshgRs72EygsNw2rHprqHEJLSwsSEhJw/fp1REVFCV0OIUPS6/WYP38+KisrIZHQe1sirLPlJmx4/Qysds8Crz8quRS7Ns1EdmK4R4+nfxVDGDNmDB5++GG89dZbQpdCCC+pqamIiIjA6dOnhS6FBLiaZivyd3o39IDO0d/GnWc8nvak4OOha+NqGhwTsaDF7ERojDE8+04xOuy+WZrQYXdh8+5ij16XKfh4mDdvHlpbW3Hu3DmhSyGEl67gozdrRCj7zlehpLoZDpdvfgYdLoYLlc3YX1zl9mMp+HiQSCTIz8+nJhciGhkZGZDL5Th//rzQpZAAxBjD9qN6jxtZ+LLandh2RO/2GzxqbuHJYDAgKysLlZWVCAoKErocQob005/+FFKpFC+++KLQpZAAc7bchI1vnBky+BoP/Tfaqy7D2dIATiqHYnwqNAuegCIqkfe11Aop3nzCvUYXGvHxpNPpMH36dOzfv1/oUgjhhZY1EKG8W2Tk1dDSeuETSJTBCJ5yHzilGu2lRbi15wUwh433tax2J/YUGd2qj4LPDV1NLoSIwYwZM9Da2orLly8LXQoJMIXlZvCZS4zJ34HYDX9GxAM/RMxjvwUAOG83wtZQwftajHVezx0UfG7Izc1FYWEhjEb33l0QIgSO42jUR0ac3emC0WTh9bXKmIndnzOXo/MTTgJpiHvr84xmi1sbW1PwuUGlUmH16tUoKCgQuhRCeKFlDWSkNbR2QCZ1L1pcNisaP9wBABgzcyVkbgafTMKhoZX/cVwUfG7atGkT3njjDbhcvj02gxBvmDNnDqqrq1FaWip0KSRA2J0M7uwj7bQ0o+7tn6Oj6gpC7rofY+c/4fY1OY6D3cm/T5OCz03Z2dlQqVT44osvhC6FkCFJpVKsXLmSpjvJiJFLObh4rt1zNN9C7Vv/DlvtdYyZvQoRD/zAo9MXGGOQS/k/joLPTRzHUZMLERWa7iS+0tLSgq+++gp/+9vfsGXLFixZsgTZ6WmwdvDryqzdtRUOUxWkY6LA7B0wHX0VpqOvoqP6mlt1OFwMkSFK3l9PB3Z54PHHH0daWhpu376N0NBQocshZFDz58+HXq9HVVUV4uLihC6HiFBXd/ClS5d6fTQ2NmLy5MlIT09Heno6cnJykJ6ejifeK0NZw9ANLs7WzlMWnC31uF34fvefK8YlQzk+jXd9Wo0acjfuK1LweSA6Ohrz58/Hnj178OSTTwpdDiGDUigUWL58Ofbt24fvf//7QpdD/JjFYsHVq1dRUlLSK+Dq6uqQlpbWHXDPPPMM0tPTkZSU1O8JIDMSm1HeaBlySYPupx8Mu2aOA7ITNe49hnZu8cz777+Pl156CSdPnhS6FEKG9P7772P79u04fvy40KUQP9De3o5r1671CbiqqiqkpKR0B1zXx4QJEyCVSnk/P9+dW7zBk51bKPg8ZLfbodVqceLECaSl8R+SEyIEq9WK2NhYOlcywNhsNuj1+j4BZzAYMGHChD4BN3HiRMjl8mFflzGG+/54HEbz8E9LH0pCuBonts53qymGgm8Ytm7dCrlcjt/97ndCl0LIkNasWYPFixfjqaeeEroU4mV2ux03btzoE3ClpaVITEzsE3CpqalQKBQ+rWnvuUr88kCJT0d9KrkUv83NQO70eLceR8E3DJcuXcLixYtRUVEBmYxulxL/tnv3bhQUFODQoUNCl0I85HQ6cfPmTVy6dKlXyN24cQPx8fF9Ai4tLU2wTfUZY1j72mkUGcw+OZpIJuGQpdPgnX+7x+0lEBR8wzRr1iy88MILWLZsmdClEDKo27dvIz4+HgaDAWPHjhW6HDIIl8uFsrKy7mDrCjm9Xo/o6GhkZGT0CrhJkyZBrVYLXXYfNc1W5Gw/gbYO74/6gpVSHNsyHzFh7gc7Bd8w/eUvf8GxY8fw7rvvCl0KIUNasWIFVq9ejXXr1gldCkFnwFVUVPQJuKtXryIiIqJPwE2ePBkhISFCl+2WwnIT1r9+htdpDXyp5FLs2uReQ0tPFHzD1NTUhMTERNy4cQORkZFCl0PIoAoKCnDgwAHayWWEMcZQWVnZJ+CuXLmCMWPG9Am4KVOmYMyYMUKX7TWF5SZs3HkGHXbXsKY9ZRIOSrkEBfmehx5AwecV69atw6xZs/DDH/5Q6FIIGZTJZEJSUhKqq6sRHBwMu9OFhtYO2J2dWz5FhijdWghMemOMoaampk/AXb58GSqVqt+A02jcW4MmVjXNVmzeXYwLlc0ejf7UCiky48KwY800xIaphlULBZ8XHDt2DD/60Y9QXFwsdCmEDIoxhvvynkDsvbloQBiMJgtkUgk4rvNcM4fTBW24GtmJGqzK0iJbp/Fo78TRjjGGW7du9Qm4S5cuQSaT9Qm49PR0RERECF224Bhj2F9chW1H9Ghss8Fqdw66yJ3jOqc1I4IVeG5xKlZOi/PKzyMFnxe4XC4kJydj3759mD59utDlENIHYwz7zldh+1E9apvaYHcB4AYe2fV8wdmSk4rc6d55wRGjhoaGXksEukLO5XL1G3Djxo0TumS/xxhDkcGMPUVGFJabYTRbIJNw4DgOjDE4XAxaTecbsNVZWmR5+Q0YBZ+XvPDCCzCbzfjv//5voUshpJeuKaaLVc0eralSyaWYGu+dKSZ/Zjab+w249vb2fgMuJiYmYN8MeNtIT7lT8HlJWVkZZsyYgaqqKiiV/HcJJ8SXzpabkO9HTQX+oKWlpc9myyUlJbh9+3afcEtPT0dcXOCOdkcrCj4vWrhwIZ555hmsWrVK6FIIwdlyEzb4WRv5SOrvRIGSkhKYTCZMmTKlT8AlJCRQwAUICj4veuutt/D3v/8dH330kdClkADn64XDR7fM85tpT4vFgitXrvQJuFu3biEtLa3PNGViYmK/JwqQwEHB50UWiwXx8fG4ePEinXtGBMMYw5pXT+Nchf9tFTUc7e3tuHr1ap9pyq4TBe4MuOTkZLdOFCCBg4LPy55++mkkJSXhZz/7mdClkADlz5sD82Gz2XDt2rU+AVdRUYHk5OQ+ATdx4kTaK5e4hYLPy06fPo3169dDr9fT/QIy4kbyOBitRoXPf7zA459zu92O69ev9wm4rhMF7gy4lJQUn58oQAIDvU3yslmzZkEmk+HUqVO49957hS6HBJhCgxmNbTZeX9ty9gBaLxyBvaECYC6EzXkMY+fy38Ozsc2GIoN5yEYXh8PRfaJAz4+uEwW6Ai4vLw/PP/880tLSqDOa+BQFn5dxHIdNmzbh9ddfp+AjI+7dIiPvLk5b7Q1IgkIgDY2Es+WW29ey2p3YU2TsDr6uEwXuPBNOr9cjJiYG6enpyMjIwPLly/GTn/wEkyZNgkrlHw0yJLDQVKcP1NbWYvLkyTAajaLbSZ2I28I/f4bShja3HnPrvf+C9fppt0d8ADBW0o6pxgPdJwpERkZ2B1zPEwWCg4Pdel5CfIlGfD4QExODe++9F//85z+Rn58vdDkkQNidLhhNlhG9ZrNTgbnz5uN73/sepkyZgtDQ0BG9PiGeoMUsPtI13UnISGlo7YBshE9WCFLIsPzRxzBr1iwKPSIaFHw+8uCDD+Lq1au4fv260KWQAGF3Mox0IzHHcbA76W4JERcKPh9RKBR4/PHHsXPnTqFLIQHA4XCgtqoSDofv1u71h7HOTYUJERNqbvGhixcv4oEHHoDBYKAdJMiwMMZQX1+PsrIylJaW9vm1uroaMbHjwa15GZDw+1m7/c3H6DBeRnvFBThb6iEflwTFuGSoU++BOnU2r+eQSzlc/s1SOryWiAo1t/hQZmYmYmJicPToUdx///1Cl0P8XFtbG8rLy/sNtrKyMiiVSiQlJSE5ORlJSUmYMWMGVq9ejeTkZCQkJEChULjV1dlhvIy2kmPdv7ffKoP9VhlkYeN4B59Wo6bQI6JDwedjXU0uFHzE6XSisrJywGBrbm5GYmJir3CbN28ekpKSkJSUhLCwsCGvkZ2oQVlj26CnWneJXL4Fkcu3ePz9cFzn9QgRG5rq9DGz2YzExESUlZUhPNz/j3IhnmOMobGxccDpyMrKSowbN65XsPX8NSYmZtinBpwtN2HjG2d8uk9nF7VCijefEMcRRYT0RCM+H9NoNFi2bBnefvttfP/73xe6HDJMVqt10OlIqVTaK9CmT5+OvLw8JCcnQ6fT+XwrrmydBhHBClhsvt+rMzJEiSwdjfiI+NCIbwR88skn+NnPfoaioqLuP7M7XWho7YDd2dkVFxmipHslfsDpdKK6unrAYDOZTNDpdP2O2pKSkqDRCB8EYj+dgRBfo+AbAU6nE4mJiXjpjfdwvlmJwnIzjCYLZFIJOA5gDHA4XdCGq5GdqMGqLC2ydRo63cFHzGbzgMFWUVGBiIiIAacjx48f7/eHmDLGsPa10ygyjK7z+AjxFgo+H2OMYd/5Kvzq3a9hZTK4JLJBGw84rvPddESwAltyUpE7PY5eXNzU3t4Og8EwYLgxxnqN0noGm06nGxUbJ/v6BPZjW+YjJizI689NyEig4POhmmYrNu8uxsWqZo+mnVRyKabGh2HHmmmIDRP/i7G3uFwu1NTUDBhs9fX1SEhIGHDUptEExmi6sNyE9a+f4X1aAx8quRS7NlFDCxE3Cj4fOVtuQv7OM+iwu4Y13SSTcFDKJSjID6wXm+bm5gGDzWAwYOzYsQMGW1xcHG0Y8C+F5SZspJ9DQnqh4POBs+UmbKB32oOy2Wz9Tkd2fW632/sNtaSkJCQmJkKtVgv9LYhG18zDhcpmj34m1QopMuNo5oGMHhR8XubreytHt8wTxYsPYwy1tbUDjtrq6uoQHx8/4KgtIiIiIKYjRwpjDPuLq7DtiB6NbTZY7U7e95qfW5yKldPoXjMZPSj4vIgxhjWvnsa5isDopmtpaek1Suv5a3l5OUJDQwcMtvj4eMhktIx0pDHGUGQwY0+RsbO72GyBTMKB4zgwxuBwMWg1nd3Fq7O0yKLuYjIKUfB50WhbP2W321FRUdFnGrLrV6vVOuh0JJ0+7/9oPSkJRBR8XsIYw31/PA6j2fc7Zmg1Knz+4wXDfifOGMOtW7cGnI6sqanB+PHjBxy1RUVF0WiAECI6NNfkJYUGMxrbbLy+ljlsMH/6OtqufgFms0IRPQGaRU9BOT6N1+Mb22woMph5Nbq0trYOOB1ZVlYGtVrdK9DuuecePPbYY0hOToZWq4VcLudVEyGEiAUFn5e8W2Tk3TFnOvoqWosPQx6lg1x3FyxXvkDdO79E3P/6G6TqoXfgt9qd2FNkRHZiOBwOB4xG44AbI7e2tvZZqL1w4cLuPwsNDR3ut04IIaJCweclheVmXkfBONua0HrhKMBJEL32RUiDx6JBIkXbpeO4XfQBxs5dN+RzMAbs++IC/vnjXFRVVSEmJqZXsC1fvrz799HR0TQdSQghPVDweYHd6YLRZOH3tQ0VgMsBaVg0pMFjAQCKmIlou3QctltlvK/pUofj0OGPkZyog0Kh8KRsQggJSNS+5QUNrR2Q8eyEc7aZAQASxf/sc8j96/Ou/8aHQiZBWHQ8hR4hhLiJgs8L7E4GvrOJ0uDOY2tctvbuP2P/+rzrv/HBcRzsTmrIJYQQd1HweYFcyvG6vwcA8kgtIJHB2VLfPcLrqNEDABTjknhfk7HOdVeEEELcQ/f4vCAyRAmH08Xra6XBGoRkLkLrNx+j7h+/gDxKB8uVk+AUKoRmLed9TYeLITLEt6d5E0LIaETB5wVyqQTacDVKG9p4fb0m52lAKoPlyhewm2ugjEuDZuGTvJYydNFq1LTDBiGEeICCz0uyEzUoa2zjNeUpkSsRseQZRCx5xqNrcVzn9QghhLiPhgxesipLC5V8ZM6AU8mlWJ2lHZFrEULIaEPB5yXZOg0igkdmaUFkiBJZOhrxEUKIJyj4vITjOGzJSYVa4dtRn0ouxZacFNqNhRBCPETB50W50+OQGRcGmcQ3oSSTcJgaH4aV0+J88vyEEBIIKPi8iOM47FgzDUq5b/5alXIJXl4znUZ7hBAyDBR8XhYbpkJB/kyvN7qo5FIU5M9ETFjQ0F9MCCFkQHQQrY8UlpuwcecZdNhdcLg8/yuWSTgo5RIU5M/kdf4eIYSQwVHw+VBNsxWbdxfjQmUz77P6elIrpMiMC8OONdMQG6byQYWEEBJ4KPh8jDGG/cVV2HZEj8Y2G6x256CL3Dmuc1ozIliB5xanYuW0OLqnRwghXkTBRwghJKBQcwshhJCAQsFHCCEkoFDwEUIICSgUfIQQQgIKBR8hhJCA8v8BZZMreDkH/CgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ism.draw_topology()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By sampling the final state of annealing for many times, we can estimate the $P$ parameter for a fixed inverse temperature $\\beta$.  By varying $\\beta$ and repeating this procedure, we can get a relation between $P$ and $\\beta$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "shots = 100  # Time of measurements for each beta.\n",
    "beta_list = np.linspace(-2, 2, 20)  # Different beta values.\n",
    "P_list = []\n",
    "\n",
    "for ism.beta in beta_list:\n",
    "    P = 0\n",
    "    for i in range(shots):\n",
    "        config = np.array(ism.measure()).reshape(3, 3)  # Reshape the qubit result to an array.\n",
    "        tmp = 0\n",
    "        for idx in range(3):\n",
    "            for jdx in range(3):\n",
    "                tmp += (-1)**(idx+jdx) * config[idx, jdx]\n",
    "        P += abs(tmp) / 9.  # P for a certain outcome.\n",
    "    P_list.append(P / shots)  # P for a certain beta."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The relation can be show in a figure.  We omit the error bars in analysis to simplify the calculation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEMCAYAAADAqxFbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAkqElEQVR4nO3deXhU5d3/8fd3shLCEgiEPSxGEEXBIOAOdSnVVq2lVdvyaFu1tfVqn6er/bVP9z6t2u3nT6u2dtFeVmqtrTwWC9oG9yiEKggIhiVIBJEQkBAgycz398cMmGACmWRmziTzeV3XXMxyn5lPTsJ8z32fc+5j7o6IiMghoaADiIhIelFhEBGRNlQYRESkDRUGERFpQ4VBRETayA46QHcVFxf72LFju7z8vn376Nu3b+ICJYhyxUe54qNc8emNuaqqqna6+5B2X3T3Hn0rLy/37qioqOjW8smiXPFRrvgoV3x6Yy5guXfwvaqhJBERaUOFQURE2lBhEBGRNlQYRESkDRUGERFpI6WFwcx+a2Y7zOyVDl43M7vNzKrNbKWZnZrKfCIikvoew++BuUd5/X1AWex2PXBnCjKJiMStqqaeOyqqqaqpD+zzH93QlJTPT+kJbu7+lJmNPUqTS4H7YsfYVprZQDMb7u7bUpNQRDJFVU09lRvrmDV+MOWlRe96PXpMP4RCxr5mZ9nmXUwe3p++edn8efnrfP3hVYQjTlbImH96KccNLSQvO4u87FD0lpPFqWMG0i8/h92NTdTta6J0UAHZWSEOtoRZUVPPss27mDamiEnD+nOwJczBlggHmyOH75eXFpGTFWL1G3tY/+ZePjhtFAB3P7mBWxevIxxxHt1cyf3Xzmr3Z+gq8xRfjyFWGB5195Paee1R4Mfu/kzs8T+Br7n78iPaXU+0R0FJSUn5ggULupynoaGBwsLCLi+fLMoVH+WKT6bliriz+6BTkG3kZxvLtjdz18tNhB1CBmUDQ4QMGlugsdlpbHH2t8Dnp+UxdWg2z9c0cPda49un5zNuQBa3rTjAih3hY37ud8/Ip7R/Fv/a0sx9a5r4v3MKGJBn3LPqAM/UHnv5295TQP9c46H1TSza1Mxv3xs9y/k7z+1n89sRIDrsc3lZDu+fkBvXOpkzZ06Vu09v77UeOSWGu/8K+BXA9OnTffbs2V1+r6VLl9Kd5ZNFueKjXPHpDbmWb97Fv17dwYQhhRT3y2Pn3oO81XCQt/a+c7tq5hguOWUEr25/m0/+4mnu+OipzD55OI899DJh3wpAxKGuOZuxg/syvE8O/fOz6d8nh/75OVw8dQRlJf14e0kFfzjjJE4ZPZD++Tnkj9nJNb9bRnNLhJysEHfPL2fS8P5ttvYPtkQ4YXg/CnKzGbtzH+Un72buScPIy85i8a5VWO0WHDDgvMklXDi5JNbbyCIvJ9rrKC8tIi87iynTD/LVg2FGD+qDmXHzmDqu/t2LNDVHyM0JcdX5pyW0x5BuhaEWGN3q8ajYcyIihy1auY3P/nFFu6/1ycliSL88hvTL49CIyOiiAv7ng1OYMnIAAB8qH8UjL71BczhCTnaIu+dPP+oXa/9c4+yyd6YVmjW+mPuvnXXUoajWxhb3ZWzxO3MazSsfxV//vTVaWLJD3HDuhKO+x+DCPAa36kjNHD+Y+6+dxQNPLEt4UYD0KwwLgRvNbAEwE9ij/QsiArDmjbfZsmsfc08azsadDYefDxlcedoYPn3ueIoL8+ib9+6vtb552Xx05pjDj2eMG8z913X+i7095aVFXf5CLi8tiquwdPQeeyfkJrwoQIoLg5k9AMwGis1sK/BtIAfA3e8CFgEXAdVAI/CJVOYTkfQSiTihkAHwiyfWs/qNt7lw8jBOn1BMfk714S3uD5WPonRwfLOMdueLPRGC/vyjSfVRSVcd43UHPpeiOCKSpvY0NvOn5Vv4Q2UN935iBuOHFPLf759Mv/xsQiFLyBa3dCzdhpJEJIPVNkT4xl9X8fCKWvY3h5k5bhCNTdGjd0YPKmjTNp23uHs6FQYRCVQk4lSs28Hvn9vM06/tJzd7K5dNHcE1Z4xj8oj+QcfLSCoMIhKIqpp67n5yAyu37mb72wcZ1j+fD5Xl8I0rZzOob3zH5EtiqTCISMpV1dTzsXsqOdAcwQy+eMHx3DB7As8+/ZSKQhrQ7KoiklLN4QiLV2+nqeWdM3ezQkZOlr6O0oV+EyKSUrcuXsf9lTXkZIXIMsjJDjFr/OCgY0krGkoSkZSaP6uUEQPymTJqoA43TVMqDCKSEpt37qN0cAGjBxVwzZnjAFQQ0pSGkkQk6Ta81cAHbn+GWxevCzqKdIIKg4gk1e7GJq69dzm5WaE28xVJ+tJQkogkTXM4wuf+uILa+v388bqZjCoqOPZCEjgVBhFJmu8/uoZnq+v4yYdPYfrYQUHHkU7SUJKIJMUfKmu47/kaPn3OeOaVjwo6jsRBhUFEEu7Z6p18Z+Fqzps0lK/OnRR0HImTCoOIJNSmnfv47P0rmDCkL7+4cipZsespSM+hwiAiCdUnJ4upowfym6tPo19+TtBxpAu081lEEqIlHMHMGDYgn3s/OSPoONIN6jGISEL84O9rue6+5TSHI0FHkW5Sj0FEEmLC0ELyskOaJbUXUGEQkW5paomQmx1i/qzSoKNIgqi0i0iXbd65jzk/WcpT698KOookkHoMItIlT7/2Fl9Y8BJNLWFKB2uqi95EPQYRiVtVTT1X//ZFdu1roqnF2dnQFHQkSSAVBhGJ25LV24l49H44EqFyY12wgSShVBhEJG579jcDENKlOXsl7WMQkbi4O89tqGPKyP7MPWm4Ls3ZC6kwiEhcXty0iy27GvnZR07h8lM1a2pvpKEkEYnLQ1VbKczLZu5Jw4KOIkmiwiAicYk4XDZtBAW5GnDorfSbFZG4/PQjp+DuQceQJFKPQUQ6bfueAwCY6RoLvZkKg4h0ypa6Rk7/8T/58/LXg44iSabCICKd0i8/m6++dxJnlRUHHUWSTPsYRKRTivrmcsPsCUHHkBRQj0FEjumV2j088lItTS26CE8mUGEQkWP67TOb+ObfXiGio5EyQkoLg5nNNbN1ZlZtZje18/oYM6sws3+b2UozuyiV+UTk3fYeaGbRK9v4wCkjyM/JCjqOpEDKCoOZZQF3AO8DJgNXmdnkI5p9E3jQ3acBVwK/TFU+EWnfolXbONAc4cPlmv4iU6SyxzADqHb3je7eBCwALj2ijQP9Y/cHAG+kMJ+ItOPPy7cyYUhfpo4eGHQUSRFL1RmMZjYPmOvu18YezwdmuvuNrdoMB5YARUBf4Hx3r2rnva4HrgcoKSkpX7BgQZdzNTQ0UFhY2OXlk0W54qNc8elsru37Itz09H4+cnwOF43PTZtcqdYbc82ZM6fK3ae3+6K7p+QGzAPuafV4PnD7EW2+CHwpdv90YA0QOtr7lpeXe3dUVFR0a/lkUa74KFd8Opvrln+s9XE3Perb9+xPbqCYnr6+Uq07uYDl3sH3aiqHkmqB0a0ej4o919qngAcB3P15IB/Q2TQiAQhHnIdX1HLu8UMo6Z8fdBxJoVQWhmVAmZmNM7NcojuXFx7RZgtwHoCZnUC0MLyVwowiEvNs9U627TnAvPLRx24svUrKCoO7twA3AouBtUSPPlptZt8zs0tizb4EXGdmLwMPANfEujwikmIhM84uK+b8yUODjiIpltIpMdx9EbDoiOe+1er+GuDMVGYSkfadVVaseZEylM58FpF3ee3NvexpbA46hgREhUFE3uUrD63k4795IegYEhDNrioi7/L9S09i7wH1GDKVCoOIvMuUUQOCjiAB0lCSiBzWEo7wnYWrWbvt7aCjSIBUGETksKdf28nvn9tMTd2+oKNIgFQYROSwh6q2UlSQw3smlQQdRQKkwiAiAOxubOLxNW9y6dSR5GbrqyGT6bcvIgA88tIbNIUjfHi6rruQ6VQYRASIDiOdMLw/J47QEUmZToVBRHh1+9usqt2jq7QJoMIgIsBDy7eSHTIunToi6CiSBlQYRDJcczjC316q5bwThjK4MC/oOJIGdOazSIY72BLhqhljmDV+cNBRJE2oMIhkuMK8bL504cSgY0ga0VCSSAbbtS967kJzOBJ0FEkjKgwiGezRlW9w3X3L2fiWpsCQd2goSSSDXTVjDMcNKWTisH5BR5E0oh6DSIaqqqnnV09tJC8nK+gokmbUYxDJQNX1YX685HnCEScvJ8T9186ivLQo6FiSJtRjEMlAa+rCtEQcB5pbIlRurAs6kqQRFQaRDNQnxwAIGeRkh3QOg7ShoSSRDFTbECEvO8RnZ0/grLIhGkaSNlQYRDJMSzhC1fYWLpg8nC+cf3zQcSQNaShJJMNUbtzF3mZ4/8nDg44iaUqFQSTD/H3VNvKyYPbEoUFHkTSlwiCSQVrCERav3s7UIVnk6/wF6YAKg0gGqdy4i137mjhtmHYvSsdUGEQyyJB+eXxs5hhOHqLegnRMhUEkg0wc1o8ffnAKuVkWdBRJYyoMIhliw1sNrNy6G3cPOoqkORUGkQxxz9MbuepXlRxs0bUX5Oi0B0okQ3xt7iQuOWWkjkaSY1KPQSRDDCzI5fQJmhNJjk2FQSQD/Pqpjfxp2ZagY0gPkdLCYGZzzWydmVWb2U0dtPmIma0xs9Vm9sdU5hPpjZrDEe5YWs3zGzS1tnROyvYxmFkWcAdwAbAVWGZmC919Tas2ZcDXgTPdvd7MdM6+SDc9t6GO3Y3NXDRFcyNJ56SyxzADqHb3je7eBCwALj2izXXAHe5eD+DuO1KYT6RXWrRyG4V52Zxz/JCgo0gPYak6ptnM5gFz3f3a2OP5wEx3v7FVm78B64EzgSzgO+7+j3be63rgeoCSkpLyBQsWdDlXQ0MDhYWFXV4+WZQrPsrVvpaI84WKRk4uzuLTp+SnTa6OKFd8upNrzpw5Ve4+vd0X3T0lN2AecE+rx/OB249o8yjwVyAHGAe8Dgw82vuWl5d7d1RUVHRr+WRRrvgoV/uWrtvhpV971Jes3t7m+aBzdUS54tOdXMBy7+B7NZVDSbXA6FaPR8Wea20rsNDdm919E9HeQ1mK8on0On9f+Qb98rI5u6w46CjSg6SyMCwDysxsnJnlAlcCC49o8zdgNoCZFQPHAxtTmFGk12gOR1i8+k3On1yik9okLikrDO7eAtwILAbWAg+6+2oz+56ZXRJrthioM7M1QAXwFXfXMXYiXfBs9U727G/mYh2NJHGK63BVM+sHvA+YAiwBnomNVXWKuy8CFh3x3Lda3Xfgi7GbiHRDv/wcLp4ynLOP1zCSxOeYhcHMRgCXAB8ETgH+CawCvgtMMrPFRIeAlrj7/uRFFZF4lJcWUV5aFHQM6YE602N4hui+gFuAJ2NDQgA/NrOBwMXAx4DLgE8kIaOIxGlLXSNZWcbIgX2CjiI90DELg7uPj/UaclsVhUOv7Qbuj91EJE3cXvEaj72ynRX/fQE5WZoSTeJzzL8YM/s80cNIN5jZi2ZWkvxYItIdn5tzHD/7yFQVBemSzvzVfJXoMNFIovsWvp/MQCLSfaWD+3LBZG3DSdd0pjAUuftCd98OfAmYk+RMItINf6is4Yk1bwYdQ3qwzhSG8KE7sX0KOsxBJE01tUS49R+vsuiVbUFHkR6sM0clFZrZdqJnLr8I5JrZcHfXX55Imnm2eidvH2jRSW3SLZ0pDIOAqbHbNGATUGNmDcBq4BV3vyFZAUWk8/6+ahv98rM5S3MjSTd05nDV3cDS2A2A2FxHJwGnEj3pTUQC1tQSYcnq7VwwuYS8bM2NJF3XmTOff0b0zOanD01/4dEL7awAVphZPzO7Ahjk7ncmM6yIdOzQMNL7T9YwknRPZ3Y+bwa+DdSa2e/M7DIzm2BmN8Smw1gPfABYl8ScInIMj66MDSMdpyu1Sfd0ZijpNuC22PQXFwFXAT8kOonezUSnyQh3/A4ikmxNLRGWrNnOhZOHkZutk9qkezo9u2psX8MfYzcRSSPPVL/FXg0jSYJo00KkFzAzzpgwmDOP09FI0n1xXY9BRNLTnIlDmTNxaNAxpJdQj0Gkh9u+5wANB1uO3VCkk1QYRHq4Wxa/ypyfLCUS6fTFFEWOSkNJIj3c/FmlnFM2hFDIgo4ivYQKg0gPN21MEdPGaG5LSRwNJYn0YA+v2MqKLfVBx5BeRoVBpIc62BLm24+s5oEXtgQdRXoZFQaRHur3z25m78EWJg7rF3QU6WVUGER6oKqaem75R3R6sp8sWUdVjYaTJHFUGER6oH+ufZNwdLJjmlsiVG6sCziR9CYqDCI9UE3dPgBCBjnZIWaNHxxwIulNdLiqSA+zbc9+Hl+7g/MmDeXU0iJmjR9MeakOV5XEUWEQ6WFu/1c17s53Lz2RUUUFQceRXkhDSSI9yOu7GvnTste58rQxKgqSNOoxiPQg/fvk8JlzJ/DxWaVBR5FeTIVBpAcZ0CeHL793YtAxpJfTUJJID/Gzx9dTsW5H0DEkA6gwiPQA+5vC/O/Lb7Bs066go0gG0FCSSA/QJzeLJf91Di1hXXNBkk89BpE0t33PAfY3hcnJCtEnNyvoOJIBVBhE0tzXH17JpXc8oyu0ScqktDCY2VwzW2dm1WZ201HafcjM3MympzKfSLqpqtlFxbq3+OC0UbpCm6RMygqDmWUBdwDvAyYDV5nZ5Hba9QO+ALyQqmwi6eqnS9ZTXJjL1WfovAVJnVT2GGYA1e6+0d2bgAXApe20+z5wM3AghdlE0s5zG3by3IY6bph9HAW5Ok5EUsfcUzNuaWbzgLnufm3s8Xxgprvf2KrNqcA33P1DZrYU+LK7L2/nva4HrgcoKSkpX7BgQZdzNTQ0UFhY2OXlk0W54tPbcrk7//PCAXbud24+pw+5WYkdRupt6yvZemOuOXPmVLl7+8P17p6SGzAPuKfV4/nA7a0eh4ClwNjY46XA9GO9b3l5uXdHRUVFt5ZPFuWKT2/LVfHqm176tUf9D89vTmygQ+/fy9ZXsvXGXMBy7+B7NZVDSbXA6FaPR8WeO6QfcBKw1Mw2A7OAhdoBLZnG3fnpkvWMKurDR6aPPvYCIgmWysKwDCgzs3FmlgtcCSw89KK773H3Yncf6+5jgUrgEm9nKEmkN1uy5k1W1e7h8+eVkZutI8ol9VL2V+fuLcCNwGJgLfCgu682s++Z2SWpyiGS7ooKcrl4ynAunzYy6CiSoVJ6qIO7LwIWHfHctzpoOzsVmUTSzYxxg5gxblDQMSSDqZ8qkiZawhF+ubSa+n1NQUeRDKfCIJImVmzZza2L11G5sS7oKJLhdNaMSJqYMW4Q//ziuYwr7ht0FMlw6jGIpIGGgy0AjB9SiJnmRJJgqTCIBOxAc5j3/vwpfvHE+qCjiAAqDCKBW/DiFmp37+e0sToSSdKDCoNIgPY3hbm9YgOzxg/ijAmDg44jAqgwiATqvuc3s7PhIF+6cKL2LUjaUGEQCcgz1Tv5+RPrmTp6oIaRJK2oMIgEoKqmnmt++yIHmiOs2fY2VTX1QUcSOUyFQSQAj7xUS0vsGs7hcEQntUla0QluIim2Y+8B/vflNzAgZJCTHWLWeO14lvShwiCSYoMKcvnw9NFMGtaPbXsOMGv8YMpLi4KOJXKYCoNIioQjzu7GJgYX5vF/Ljoh6DgiHdI+BpEUueUfr/KB//eMZk+VtKceg0iKXDJ1BIV52RT1zQ06ishRqTCIJNkbu/czYmAfThwxgBNHDAg6jsgxaShJJIle3xvhwp8/xT1Pbww6ikinqTCIJMmOvQf4edUB+uZlcfHJw4OOI9JpGkoSSYL9TWGuu3c5Dc3Ow9edxvABfYKOJNJp6jGIJFgk4nzxwZdYWbuHG07J46SR2q8gPYt6DCIJdsvidTz2yna+efEJHBfeEnQckbipxyCSQAte3MJdT27gYzPH8KmzxgUdR6RLVBhEEuTZ6p1882+vcHZZMd+95ERdX0F6LBUGkQQp6Z/P7IlDueNjp5Kdpf9a0nNpH4NINzU2tdAnJ4vjhhZyz9XTg44j0m3arBHphoMtYeb/5kW+s3B10FFEEkY9BpEuqqqp5/kNOxlfXMBMXU9BehEVBpEuqKqp56O/rqQ5HCE3O8SVM0qDjiSSMBpKEumCnz++noMtESIOzS26NKf0LuoxiMTpF0+s55nqnWSZAa5Lc0qvo8Ig0knuzs8fX89t/6pmXvkorjhtNC9u2qVLc0qvo8Ig0gnuzk+WrOOOig1cMX00P7p8CqGQcdrYQUFHE0k4FQaRY3B3bv7HOu56cgNXzRjDDy87iVBIZzVL76XCIHIU7s7/LFrLr5/exMdnjeF7l6goSO+X0qOSzGyuma0zs2ozu6md179oZmvMbKWZ/dPMdAygBOrZ6jp+/fQmrj69lO9fqqIgmSFlPQYzywLuAC4AtgLLzGyhu69p1ezfwHR3bzSzG4BbgCtSlVHkSGeVFXPvJ2dwTlmxJsWTjJHKHsMMoNrdN7p7E7AAuLR1A3evcPfG2MNKYFQK84kA0Qvt/PixV3mldg8A5x4/REVBMoq5e2o+yGweMNfdr409ng/MdPcbO2h/O7Dd3X/QzmvXA9cDlJSUlC9YsKDLuRoaGigsLOzy8smiXPFJZK63m5zvPrefs0Zm88Gy3LTJlUjKFZ/emGvOnDlV7t7+rI/unpIbMA+4p9Xj+cDtHbT9ONEeQ96x3re8vNy7o6KiolvLJ4tyxScRucLhiIfDEXd339Vw0CORSLffszevr2RQrvh0Jxew3Dv4Xk3lUUm1wOhWj0fFnmvDzM4HvgGc6+4HU5RNMlw44tz0l5VkhYwfXT6For7d6ymI9GSp3MewDCgzs3FmlgtcCSxs3cDMpgF3A5e4+44UZpMMFo44X3noZf5ctZVhA/K1P0EyXsoKg7u3ADcCi4G1wIPuvtrMvmdml8Sa3QoUAn82s5fMbGEHbyeSEOGI8+U/v8zDK2r54gXH85/nHx90JJHApfQEN3dfBCw64rlvtbp/firzSGZ7cVMd335kNWu37+Ur753I5+YcF3QkkbSgM58lo7g7W3Y18uCy17nzyQ1EHLJDptlRRVpRYZBezd3ZtHMfDkwYUsjmukbm/GTpu9pUbqzTDKkiMbpQj/Qq7s5rb+7l+Q3vXDhn3l3Pc0dFNQBjBxfwo8un8IsrppKfEyLL0PUURI6gHoP0WIeuuTxiYB9W1DTz4P1VvLBxF3X7mhg7uIClX5mDmfGLK6YyZlABAGbGVTPGADB6UAGVG+t0PQWRI6gwSI/j7txfuYVvLXyFSKsT90cO3MO5xw9h5vhBzBz3Tg/gnOOHtPs+5aVFKggi7VBhkB7n9n9V89PH1x9+bMDcsdnc+Zn3BBdKpBfRPgZJe28faOZnS9bx7y31AFx08nA+fc548rOj+wjyckKUl2gbRyRR9L9J0tbuxiYGFuSSHTLuq6yhb14208YUMWFIIV+/6AQuPHHY4X0Eeze9HHRckV5DhUHSSlNLhMde2cbvnt3M3gPNPP5f51KQm81TX51D//ycNm1b7yNYuimItCK9kwqDpIWdDQd54IUt/KGyhh17DzK+uC9XnzGWsDsh7F1FQUSSR4VBuqyqpr5bh3tW1dSz8KVatuxq5NkNdTS1RDjn+CHcPG8s55YN0WU0RQKiwiCd9sxrO1lfH2Y2sHzzLj589/O4R48KKh1cQEHuu/+czp04hK/NnQTA5b98loumDOfas8fzbPVOPn7PCxw62nTuicP48nsnctzQ9LsYikimyejCUFVTz6Mbmug3rr7LW7zd3WJOh+VPHjWA4sI8ttbvZ2t9I7X1+6P3dzcyrH8+91x9GgA/+Psa+kSauR54YdMuDl38z4GsUIgRA/u86zMGFbxzXYPhA/vQv090SGj55l2Hi0KWwZRRA1QURNJExhaGqpp6rvpVJU3hCH957TlGFvWhT05Wu21//KEplJcO4sn1b/GDR9fw6/+YTt2+Jq781fM0hx2Doy7/6/+Yztjivjy8Yit3Lt3AXz57Bq+92cAVdz9PS6T95fc1NtJ3xZMA/OWzZ9A/P4d7nt7IwytqWfSFs6mqqT/q8ofkZIVY9IWzAbj5H9HrGP/hUzOpqqnnI3c9T7idS7vm54QYVVTAyIF9mDSs/+Hnb//oqaz59zIAZo0fTH5OiOaWCDnZIW6Zd/Ixi9MdHz318P2zyoZw55MbDi+vKSlE0kfGFobKjXW0RCJAdIu3IDerwy3WPjnR1VSYl01ZSSG52aHo8mHv1PK52dHTRQYW5FBWUkiWGZUb6whHOl5+x479DI09zopdOKa4MO9wm2Mtf0h26J1TVUr65dEwuO/h5SOxomDAB04ZwafOGseooj4M6pvb7sVqjhtayNa86PPlpUXcf+2sLvdYuru8iCRPxhaGWeMHk5sdoqk5Qm5OiB9dfuwt3ujhkeWHl89rtcXcmeXfM6mE90wq6dTyS5cuZfbs8jbLXzZtJJdNG9nlz7/mzHFtfv7Wy199xlhOGT3wqMsfqbtTSmhKCpH0lLGF4dAW6wNPLOOq809L+RZvT19eRHqvjC0MEP1y3Dsht8tfikFvMQe9vIj0TporSURE2lBhEBGRNlQYRESkDRUGERFpQ4VBRETaUGEQEZE2zNuZEqEnMbO3gJpuvEUxsDNBcRJJueKjXPFRrvj0xlyl7t7uBdF7fGHoLjNb7u7Tg85xJOWKj3LFR7nik2m5NJQkIiJtqDCIiEgbKgzwq6ADdEC54qNc8VGu+GRUrozfxyAiIm2pxyAiIm2oMIiISBsZVxjM7FYze9XMVprZX81sYAft5prZOjOrNrObUpDrw2a22swiZtbh4WdmttnMVpnZS2a2PI1ypXp9DTKzx83stdi/7c4fbmbh2Lp6ycwWJjHPUX9+M8szsz/FXn/BzMYmK0ucua4xs7daraNrU5Dpt2a2w8xe6eB1M7PbYplXmtmp7bULINdsM9vTal19K0W5RptZhZmtif1f/EI7bRK7ztw9o27AhUB27P7NwM3ttMkCNgDjgVzgZWByknOdAEwElgLTj9JuM1CcwvV1zFwBra9bgJti929q7/cYe60hBevomD8/8Fngrtj9K4E/pUmua4DbU/X3FPvMc4BTgVc6eP0i4DGiV52dBbyQJrlmA4+mcl3FPnc4cGrsfj9gfTu/x4Sus4zrMbj7EndviT2sBEa102wGUO3uG929CVgAXJrkXGvdfV0yP6MrOpkr5esr9v73xu7fC1yW5M87ms78/K3zPgScZ+1dWDv1uVLO3Z8Cdh2lyaXAfR5VCQw0s+FpkCsQ7r7N3VfE7u8F1gIjj2iW0HWWcYXhCJ8kWmWPNBJ4vdXjrbz7FxEUB5aYWZWZXR90mJgg1leJu2+L3d8OlHTQLt/MlptZpZldlqQsnfn5D7eJbZjsAQYnKU88uQA+FBt+eMjMRic5U2ek8/+/083sZTN7zMxOTPWHx4YgpwEvHPFSQtdZr7y0p5k9AQxr56VvuPsjsTbfAFqA+9MpVyec5e61ZjYUeNzMXo1t6QSdK+GOlqv1A3d3M+vouOvS2PoaD/zLzFa5+4ZEZ+3B/hd4wN0PmtmnifZq3hNwpnS1gujfU4OZXQT8DShL1YebWSHwF+A/3f3tZH5WrywM7n7+0V43s2uA9wPneWyA7gi1QOstp1Gx55Kaq5PvURv7d4eZ/ZXocEG3CkMCcqV8fZnZm2Y23N23xbrMOzp4j0Pra6OZLSW6tZXowtCZn/9Qm61mlg0MAOoSnCPuXO7eOsM9RPfdBC0pf0/d1frL2N0XmdkvzazY3ZM+uZ6Z5RAtCve7+8PtNEnoOsu4oSQzmwt8FbjE3Rs7aLYMKDOzcWaWS3RnYdKOaOksM+trZv0O3Se6I73dIyhSLIj1tRC4Onb/auBdPRszKzKzvNj9YuBMYE0SsnTm52+ddx7wrw42SlKa64hx6EuIjl8HbSHwH7EjbWYBe1oNGwbGzIYd2i9kZjOIfn8mu7gT+8zfAGvd/WcdNEvsOkv1Hvagb0A10bG4l2K3Q0eKjAAWtWp3EdG9/xuIDqkkO9cHiY4LHgTeBBYfmYvo0SUvx26r0yVXQOtrMPBP4DXgCWBQ7PnpwD2x+2cAq2LraxXwqSTmedfPD3yP6AYIQD7w59jf34vA+GSvo07m+lHsb+lloAKYlIJMDwDbgObY39angM8An4m9bsAdscyrOMpReinOdWOrdVUJnJGiXGcR3be4stX31kXJXGeaEkNERNrIuKEkERE5OhUGERFpQ4VBRETaUGEQEZE2VBhERKQNFQYREWlDhUFERNpQYRBJgtjZxo+ZWb2Z1ZrZJ4LOJNJZKgwiyfEQ8DhQDFwHfDPYOCKd1ysn0RMJkpmdDAz22Lw2sel13go0lEgc1GMQSbwzgWfMLGRm5cDPgDsDziTSaSoMIok3FVhOdFK65UAj8NcgA4nEQ4VBJPGmEp3yeg5wHNHLRd4cZCCReKgwiCSQmWUBJwD/dveIR68W92zAsUTiosIgklgTgQLgfWaWZWZTic7rf2+gqUTioMIgkljTiF4l7qfAbuD3wOfdvTLATCJx0eGqIok1FXjA3X8UdBCRrlKPQSSxppEe100W6TIVBpHEOgV4NegQIt2haz6LiEgb6jGIiEgbKgwiItKGCoOIiLShwiAiIm2oMIiISBsqDCIi0oYKg4iItPH/ARNs0A4But8KAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "_ = plt.plot(beta_list, P_list, '.-.')\n",
    "_ = ax.set_ylabel(r'$\\langle P \\rangle$', fontsize=12)\n",
    "_ = ax.set_xlabel(r'$\\beta$', fontsize=12)\n",
    "_ = ax.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Artificial distribution\n",
    "\n",
    "Generate a given distribution is a hard problem on classical problem, since any distribution is made from the built-in pseudo random number generator.  This generator is not really random, if you are dealing a long period process, this could be a disaster.\n",
    "\n",
    "Can a quantum annealer output real random number? Yes, and it can do more than just the uniform random numbers.  We will follow a concrete case below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Problem\n",
    "\n",
    "Try to make a probabilistic distribution which returns $1,2,3$ with equal probabilities, but never returns $0$ and other numbers.  In other words, it means a distribution as \n",
    "\n",
    "$$ p(1) = p(2) = p(3) = 1/3 $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Solve it\n",
    "\n",
    "If we create an annealer with such a Hamiltonian\n",
    "\n",
    "$$ diag(3,-1,-1,-1) $$\n",
    "\n",
    "The final Gibbs state will be \n",
    "\n",
    "$$ e^{-3\\beta},e^{1\\beta},e^{1\\beta},e^{1\\beta} $$\n",
    "\n",
    "By increasing $\\beta$, the probability to return the first eigenstate will be supressed a lot.  The latter three eigenstates are $|1,-1\\rangle,|-1,1\\rangle,|-1,-1\\rangle$, which by some mapping can be transfered to $|01\\rangle, |10\\rangle, |11\\rangle$ , corresponding to $1,2,3$.\n",
    "\n",
    "Actually, this Hamiltonian is just \n",
    "\n",
    "$$ H = Z_1 + Z_2 + Z_1 Z_2 $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "beta = 10.  # Higher beta can suppress higher energy state.\n",
    "num_qubits = 2\n",
    "dis = Annealer(beta, num_qubits)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "connectin = [\n",
    "    [1, 1],\n",
    "    [0, 1]\n",
    "]\n",
    "dis.set_connection(connectin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's measure this annealer to test the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Got 0 for 0 times\n",
      "Got 1 for 312 times\n",
      "Got 2 for 354 times\n",
      "Got 3 for 334 times\n"
     ]
    }
   ],
   "source": [
    "count = [0] * 4\n",
    "for i in range(1000):\n",
    "    output = np.array(dis.measure())\n",
    "    # A proper transformation for notation is involved below.\n",
    "    output[0] = (output[0] + 1) / 2.\n",
    "    output[1] = (output[1] + 1) / 2.\n",
    "    count[int(3 - output[0] * 2 - output[-1])] += 1\n",
    "    \n",
    "print('Got 0 for '+str(count[0])+' times')\n",
    "print('Got 1 for '+str(count[1])+' times')\n",
    "print('Got 2 for '+str(count[2])+' times')\n",
    "print('Got 3 for '+str(count[3])+' times')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAD4CAYAAAD2FnFTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAATWklEQVR4nO3df7Ad5X3f8fcnAoMTMwbCLVVAtsAl48FuLMgtJnGmQ7CdYDqNyDRlIG3AHjqKJ7hjt2kbnHQcuy1TO01M7KRDRzbUouPaUP8oqovbypjW9UyBCCJjfhhb2DBIFUjxDwxxQgP59o/zCA6X50r3oHv3HHTfr5mds/vss+d870qrj/bZPXtTVUiStNAPTbsASdJsMiAkSV0GhCSpy4CQJHUZEJKkriOmXcByOeGEE2r9+vXTLkOSXlTuuOOOP6mqud66wyYg1q9fz/bt26ddhiS9qCR5aLF1DjFJkroMCElSlwEhSeoyICRJXQaEJKnLgJAkdRkQkqQuA0KS1DVIQCQ5OsntSb6S5J4k72vtH0vyrSQ72rShtSfJh5PsTHJXkjOHqFOS9Kyhvkn9JHBuVT2R5Ejgy0k+39b906r61IL+bwFOa9Prgavbq3RYWn/Ff512CVP14Pv/1rRLUMcgZxA18kRbPLJNB/pVdhuB69p2twLHJlm70nVKkp412DWIJGuS7AD2Atuq6ra26so2jHRVkqNa20nAw2Ob72ptC99zU5LtSbbv27dvJcuXpFVnsICoqqeragNwMnBWktcC7wZeDfwN4HjgNyZ8z81VNV9V83Nz3YcRSpJeoMHvYqqq7wG3AOdV1Z42jPQk8O+Bs1q33cC6sc1Obm2SpIEMdRfTXJJj2/xLgTcDX9t/XSFJgAuAu9smW4FL2t1MZwOPVdWeIWqVJI0MdRfTWmBLkjWMQumGqvpcki8mmQMC7ADe3vrfBJwP7AR+ALxtoDolSc0gAVFVdwFndNrPXaR/AZevdF2SpMX5TWpJUpcBIUnqMiAkSV1DXaSWpBXjo0pW5lElnkFIkroMCElSlwEhSeoyICRJXQaEJKnLgJAkdRkQkqQuA0KS1GVASJK6DAhJUpcBIUnqMiAkSV0+rE/LwoelrczD0qRp8gxCktRlQEiSugwISVLXIAGR5Ogktyf5SpJ7kryvtZ+S5LYkO5Ncn+Qlrf2otryzrV8/RJ2SpGcNdQbxJHBuVb0O2ACcl+Rs4APAVVX114DvApe1/pcB323tV7V+kqQBDRIQNfJEWzyyTQWcC3yqtW8BLmjzG9sybf0bk2SIWiVJI4Ndg0iyJskOYC+wDXgA+F5VPdW67AJOavMnAQ8DtPWPAT/aec9NSbYn2b5v374V/gkkaXUZLCCq6umq2gCcDJwFvHoZ3nNzVc1X1fzc3Nyhvp0kaczgdzFV1feAW4CfAo5Nsv/LeicDu9v8bmAdQFv/cuDbw1YqSavbUHcxzSU5ts2/FHgzcB+joPil1u1S4MY2v7Ut09Z/sapqiFolSSNDPWpjLbAlyRpGoXRDVX0uyb3AJ5P8K+CPgWta/2uA/5BkJ/Ad4KKB6pQkNYMERFXdBZzRaf8mo+sRC9v/HPi7A5QmSVqE36SWJHUZEJKkLgNCktRlQEiSugwISVKXASFJ6jIgJEldBoQkqcuAkCR1GRCSpC4DQpLUZUBIkroMCElSlwEhSeoyICRJXQaEJKnLgJAkdRkQkqQuA0KS1GVASJK6BgmIJOuS3JLk3iT3JHlna39vkt1JdrTp/LFt3p1kZ5L7k/z8EHVKkp51xECf8xTw61V1Z5JjgDuSbGvrrqqq3x3vnOR04CLgNcCPAV9I8uNV9fRA9UrSqjfIGURV7amqO9v848B9wEkH2GQj8MmqerKqvgXsBM5a+UolSfsNfg0iyXrgDOC21vSOJHcluTbJca3tJODhsc12ceBAkSQts0EDIsnLgE8D76qq7wNXA68CNgB7gN+b8P02JdmeZPu+ffuWu1xJWtUGC4gkRzIKh49X1WcAqurRqnq6qv4S+AjPDiPtBtaNbX5ya3uOqtpcVfNVNT83N7eyP4AkrTJD3cUU4Brgvqr64Fj72rFuvwjc3ea3AhclOSrJKcBpwO1D1CpJGhnqLqY3AL8CfDXJjtb2m8DFSTYABTwI/CpAVd2T5AbgXkZ3QF3uHUySNKxBAqKqvgyks+qmA2xzJXDlihUlSTogv0ktSeoyICRJXQaEJKnLgJAkdRkQkqQuA0KS1GVASJK6DAhJUpcBIUnqMiAkSV0GhCSpy4CQJHUZEJKkLgNCktRlQEiSugwISVKXASFJ6jIgJEldSw6IJP9kkfZ/vHzlSJJmxSRnEO9ZpP2fL0chkqTZcsTBOiQ5t82uSfKzQMZWnwo8vhKFSZKm66ABAVzTXo8Grh1rL+AR4B8e7A2SrAOuA05s222uqg8lOR64HlgPPAhcWFXfTRLgQ8D5wA+At1bVnUv5gSRJy+OgAVFVpwAkua6qLnmBn/MU8OtVdWeSY4A7kmwD3grcXFXvT3IFcAXwG8BbgNPa9Hrg6vYqSRrIkq9BjIdDkh8an5aw7Z79ZwBV9ThwH3ASsBHY0rptAS5o8xuB62rkVuDYJGuXWqsk6dBNchfTmUn+T5I/Bf6iTU+11yVLsh44A7gNOLGq9rRVjzAagoJReDw8ttmu1rbwvTYl2Z5k+759+yYpQ5J0EJPcxbQFuAWYZ3Rx+lTglPa6JEleBnwaeFdVfX98XVUVo+sTS1ZVm6tqvqrm5+bmJtlUknQQS7lIvd8rgd9q/5BPLMmRjMLh41X1mdb8aJK1VbWnDSHtbe27gXVjm5/c2iRJA5nkDOKzwM+9kA9pdyVdA9xXVR8cW7UVuLTNXwrcONZ+SUbOBh4bG4qSJA1gkjOIo4HPJvkyo+sFz1jC3U1vAH4F+GqSHa3tN4H3AzckuQx4CLiwrbuJ0S2uOxnd5vq2CeqUJC2DSQLi3jZNrKq+zHO/YDfujZ3+BVz+Qj5LkrQ8lhwQVfW+lSxEkjRblhwQY4/ceJ6q+uLylCNJmhWTDDFds2B5DngJo+8oLPlWV0nSi8MkQ0ynjC8nWcPoSa4+rE+SDkMv+BcGVdXTwJXAP1u+ciRJs+JQf6Pcm4G/XI5CJEmzZZKL1A/z3Edh/DCj70b82nIXJUmavkkuUv/9Bct/Cnx94TOVJEmHh0kuUv8vGD3qm9FTVx+tKoeXJOkwNcnjvo9Jch3wZ4wenPdnSbYkefmKVSdJmppJLlL/AfAjwF8HXtpefxj48ArUJUmaskmuQZwHnFpVP2jLX0/yNuCB5S9LkjRtk5xB/Dmjb0+POwF4cvnKkSTNiknOID4KbEvyQUaP5n4l8I+Aj6xEYZKk6ZokIK5kdHH67wE/Bvxf4HeqauEzmiRJh4FJhpg+BNxfVW+qqtOr6k3AfUl+f2VKkyRN0yQBcTGwfUHbHcAvL185kqRZMUlAFLBmQduaCd9DkvQiMck/7v8b+Jftm9T7v1H93tYuSTrMTHKR+p3A54A9SR4CXgHsAf72ShQmSZquJZ9BVNUu4ExgI/BvgAuAn2ztB5Tk2iR7k9w91vbeJLuT7GjT+WPr3p1kZ5L7k/z8JD+QJGl5THIGQXs4361tmsTHgD8ErlvQflVV/e54Q5LTgYuA1zC6nfYLSX68/YIiSdJABrnAXFVfAr6zxO4bgU9W1ZNV9S1gJ3DWihUnSeqa9h1I70hyVxuCOq61nQQ8PNZnV2t7niSbkmxPsn3fvn0rXaskrSrTDIirgVcBGxhd7P69Sd+gqjZX1XxVzc/NLXxMlCTpUEwtIKrq0ap6ul3X+AjPDiPtBtaNdT25tUmSBjS1gEiydmzxF4H9dzhtBS5KclSSU4DTgNuHrk+SVruJ7mJ6oZJ8AjgHOCHJLuC3gXOSbGD0De0HgV8FqKp7ktwA3As8BVzuHUySNLxBAqKqLu40L/oU2Kq6ktHTYyVJUzLtu5gkSTPKgJAkdRkQkqQuA0KS1GVASJK6DAhJUpcBIUnqMiAkSV0GhCSpy4CQJHUZEJKkLgNCktRlQEiSugwISVKXASFJ6jIgJEldBoQkqcuAkCR1GRCSpC4DQpLUNUhAJLk2yd4kd4+1HZ9kW5JvtNfjWnuSfDjJziR3JTlziBolSc811BnEx4DzFrRdAdxcVacBN7dlgLcAp7VpE3D1QDVKksYMEhBV9SXgOwuaNwJb2vwW4IKx9utq5Fbg2CRrh6hTkvSsaV6DOLGq9rT5R4AT2/xJwMNj/Xa1tudJsinJ9iTb9+3bt3KVStIqNBMXqauqgHoB222uqvmqmp+bm1uByiRp9ZpmQDy6f+iove5t7buBdWP9Tm5tkqQBTTMgtgKXtvlLgRvH2i9pdzOdDTw2NhQlSRrIEUN8SJJPAOcAJyTZBfw28H7ghiSXAQ8BF7buNwHnAzuBHwBvG6JGSdJzDRIQVXXxIqve2OlbwOUrW5Ek6WBm4iK1JGn2GBCSpC4DQpLUZUBIkroMCElSlwEhSeoyICRJXQaEJKnLgJAkdRkQkqQuA0KS1GVASJK6DAhJUpcBIUnqMiAkSV0GhCSpy4CQJHUZEJKkLgNCktRlQEiSuo6YdgFJHgQeB54Gnqqq+STHA9cD64EHgQur6rvTqlGSVqNZOYP42araUFXzbfkK4OaqOg24uS1LkgY0KwGx0EZgS5vfAlwwvVIkaXWahYAo4H8kuSPJptZ2YlXtafOPACdOpzRJWr2mfg0C+Jmq2p3krwDbknxtfGVVVZLqbdgCZRPAK17xipWvVJJWkamfQVTV7va6F/gscBbwaJK1AO117yLbbq6q+aqan5ubG6pkSVoVphoQSX4kyTH754GfA+4GtgKXtm6XAjdOp0JJWr2mPcR0IvDZJPtr+Y9V9d+S/BFwQ5LLgIeAC6dYoyStSlMNiKr6JvC6Tvu3gTcOX5Ekab+pX4OQJM0mA0KS1GVASJK6DAhJUpcBIUnqMiAkSV0GhCSpy4CQJHUZEJKkLgNCktRlQEiSugwISVKXASFJ6jIgJEldBoQkqcuAkCR1GRCSpC4DQpLUZUBIkroMCElSlwEhSeqa6YBIcl6S+5PsTHLFtOuRpNVkZgMiyRrg3wJvAU4HLk5y+nSrkqTVY2YDAjgL2FlV36yq/wd8Etg45ZokadU4YtoFHMBJwMNjy7uA1493SLIJ2NQWn0hy/yLvdQLwJ8te4fKZ9fpg9mucan35wEG7uP8OwP13aA5x/71ysY1mOSAOqqo2A5sP1i/J9qqaH6CkF2TW64PZr9H6Do31HZrDtb5ZHmLaDawbWz65tUmSBjDLAfFHwGlJTknyEuAiYOuUa5KkVWNmh5iq6qkk7wD+O7AGuLaq7nmBb3fQYagpm/X6YPZrtL5DY32H5rCsL1W13IVIkg4DszzEJEmaIgNCktR1WAZEkuOTbEvyjfZ63CL9nk6yo00rfgH8YI8OSXJUkuvb+tuSrF/pmias761J9o3ts38wcH3XJtmb5O5F1ifJh1v9dyU5c8bqOyfJY2P77z0D1rYuyS1J7k1yT5J3dvpMbf8tsb6p7b/2+UcnuT3JV1qN7+v0mdoxvMT6JjuGq+qwm4DfAa5o81cAH1ik3xMD1rQGeAA4FXgJ8BXg9AV9fg34d23+IuD6GavvrcAfTvHP9W8CZwJ3L7L+fODzQICzgdtmrL5zgM9Nad+tBc5s88cAX+/8+U5t/y2xvqntv/b5AV7W5o8EbgPOXtBnmsfwUuqb6Bg+LM8gGD2SY0ub3wJcML1SnrGUR4eM1/0p4I1JMkP1TVVVfQn4zgG6bASuq5FbgWOTrB2muiXVNzVVtaeq7mzzjwP3MXpawbip7b8l1jdVbb880RaPbNPCu3ymdgwvsb6JHK4BcWJV7WnzjwAnLtLv6CTbk9ya5IIVrqn36JCFB8AzfarqKeAx4EdXuK7nfXbTqw/g77Thh08lWddZP01L/Rmm6afaEMDnk7xmGgW0YY8zGP0Pc9xM7L8D1AdT3n9J1iTZAewFtlXVovtwCsfwUuqDCY7hF21AJPlCkrs703P+11uj86rFUvSVNfr6+S8Dv5/kVStd94vcfwHWV9VPANt49n9KWpo7Gf2dex3wB8B/HrqAJC8DPg28q6q+P/TnH8xB6pv6/quqp6tqA6MnO5yV5LVD13AgS6hvomP4RRsQVfWmqnptZ7oReHT/qXF73bvIe+xur98E/iej/7WslKU8OuSZPkmOAF4OfHsFa+p+dvO8+qrq21X1ZFv8KPCTA9W2VDP9eJaq+v7+IYCqugk4MskJQ31+kiMZ/eP78ar6TKfLVPffweqb9v5bUMv3gFuA8xasmuYx/IzF6pv0GH7RBsRBbAUubfOXAjcu7JDkuCRHtfkTgDcA965gTUt5dMh43b8EfLGdAQ3hoPUtGI/+BUbjxLNkK3BJuxvnbOCxsaHGqUvyV/ePRyc5i9HxN8g/Hu1zrwHuq6oPLtJtavtvKfVNc/+1z5xLcmybfynwZuBrC7pN7RheSn0TH8NDXWEfcmI05ncz8A3gC8DxrX0e+Gib/2ngq4zu1vkqcNkAdZ3P6O6MB4Dfam3/AviFNn808J+AncDtwKkD77eD1fevgXvaPrsFePXA9X0C2AP8BaPx8cuAtwNvb+vD6JdMPdD+TOdnrL53jO2/W4GfHrC2n2E01HoXsKNN58/K/ltifVPbf+3zfwL441bj3cB7WvtMHMNLrG+iY9hHbUiSug7XISZJ0iEyICRJXQaEJKnLgJAkdRkQkqQuA0KS1GVASJK6/j8VKRUJUsN0MgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "_ = plt.bar(range(4),count)\n",
    "_ = plt.ylabel('count',fontsize=12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the figure, we can confirm the desired distribution is produced on the annealer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Think\n",
    "\n",
    "An interesting problem is can _phalanx_ generate any distribution on natural numbers.  This means can _phalanx_ create a diagonal Hamiltonian holds any possible real numbers on its diagonal places.  For example, can we set an annealer with such a Hamiltonian\n",
    "\n",
    "$$ diag(4,0,0,0) $$\n",
    "\n",
    "The answer is no.  This is limited by our Hamiltonian's formalism, the form is lack of a part corresponding to identity matrix."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
